Ephedra regional gradient (ERG) analyses

ecoblender

alex.filazzola

Abstract

An under-examined component of the shrub-annual relationship is how regional drivers, such as climate, may alter the sign or magnitude of positive interactions. Interspecific interactions between plants have been shown to be strongly linked to climate, particularly temperature and precipitation. The stress-gradient hypothesis (SGH) predicts that higher abiotic stress (i.e. for deserts, increasing temperature and reduced precipitation) will increase the frequency of positive interactions among shrubs and their annual understory. Regional climate gradients also have indirect effects on plant composition, such as determining consumer abundance and soil nutrient composition. Nutrient availability is particularly affected by precipitation because of altered decomposition rates of organic matter and mineralization. Therefore, the strength of facilitation and operating mechanism of a shrub on the annual plant community may change along a regional gradient.

Hypothesis

We tested the hypothesis that positive interactions among shrubs and annual plants will increase with abiotic stress and reduce nutrient availability along a regional gradient of aridity.

Methods

load GAMs

Year<-c("2016","2016","2016","2016","2016","2016","2016","2017","2017","2017","2017","2017","2017","2017")
Site<-c("Barstow","Cuyama","HeartofMojave","PanocheHills","SheepholeValley","Tecopa","TejonRanch","Barstow","Cuyama","HeartofMojave","PanocheHills","SheepholeValley","Tecopa","TejonRanch")
gamsq <-c(7992.36,6789.76,7903.21,5715.36,7779.24,7482.25,7447.69,6872.41,5882.89,7157.16,5595.04,7293.16,7106.49,6955.56)
gams.richard <-c(89.4,82.4,88.9,75.6,88.2,86.5,86.3,82.9,76.7,84.6,74.8,85.4,84.3,83.4)
gams <- c(89.6,71,84.7,62.7,88.3,80.6,80.7,63.6,51.4,73.4,49.5,69.4,68.3,65.8)

gams.data <- data.frame(Year,Site,gamsq,gams)

Weather for growing season

Climate patterns within study

season1.sjd <- season1 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season1.mnp <- season1 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.sjd <- season2 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.mnp <- season2 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))

## Rain vs Temperature in 2016
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot1 <- barplot(height=season1.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot1[,1], season1.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot1 <- barplot(height=season1.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
points(plot1[,1], season1.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)

## Rain vs Temperature in 2017
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot2 <- barplot(height=season2.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot2[,1], season2.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature San Joaquin (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot2 <- barplot(height=season2.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
points(plot2[,1], season2.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature Mojave (°C)", 4, line=3)

### 2016 The rain was inconsistent and mostly absent in the Mojave. This resulted in low germination and producitivty at the southern sites
### 2017 The rain was more plentiful, but in the northern sites, there appears to be a frost period after the majority of the rainfall. Need to check number of frost days

season1.frost <- season1 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season1.frost)
##              Site frost.days
## 1         Barstow  24.725275
## 2          Cuyama  29.670330
## 3   HeartofMojave  25.824176
## 4    PanocheHills  10.439560
## 5 SheepholeValley   6.043956
## 6          Tecopa  23.626374
## 7      TejonRanch  50.549451
season2.frost <- season2 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season2.frost)
##              Site frost.days
## 1         Barstow  17.679558
## 2          Cuyama  19.889503
## 3   HeartofMojave  16.574586
## 4    PanocheHills  14.917127
## 5 SheepholeValley   1.785714
## 6          Tecopa  25.966851
## 7      TejonRanch  35.911602
## Both years had comparable number of frost days

## Compare number of consecutive frost days (i.e. frost periods)
season1[,"frost"] <- ifelse(season1$min.temp<0, -99,season1$min.temp) ## identified days below freezing
season2[,"frost"] <- ifelse(season2$min.temp<0, -99,season2$min.temp) ## identified days below freezing
count.consec <- function(x) {max(rle(as.character(x))$lengths)}

season1.frost <- season1 %>% group_by(Site)  %>% summarize(count.consec(frost))
data.frame(season1.frost)
##              Site count.consec.frost.
## 1         Barstow                  10
## 2          Cuyama                  10
## 3   HeartofMojave                  12
## 4    PanocheHills                   5
## 5 SheepholeValley                   7
## 6          Tecopa                  11
## 7      TejonRanch                  14
season2.frost <- season2 %>% group_by(Site) %>% summarize(count.consec(frost))
data.frame(season2.frost)
##              Site count.consec.frost.
## 1         Barstow                   5
## 2          Cuyama                   6
## 3   HeartofMojave                   7
## 4    PanocheHills                   6
## 5 SheepholeValley                   3
## 6          Tecopa                   7
## 7      TejonRanch                  13
## compare only after plants have germinated
season1.frost <- season1 %>% group_by(Site) %>% filter(year>2015) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100, avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season1.frost)
##              Site frost.days avg.min.temp
## 1         Barstow  12.396694     5.750413
## 2          Cuyama  11.570248     3.352893
## 3   HeartofMojave  16.528926     4.542149
## 4    PanocheHills   2.479339     6.777686
## 5 SheepholeValley   2.479339     9.394167
## 6          Tecopa  11.570248     6.342149
## 7      TejonRanch  33.884298     1.438017
season2.frost <- season2 %>% group_by(Site) %>%  filter(year>2016) %>%  summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100,avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season2.frost)
##              Site frost.days avg.min.temp
## 1         Barstow  11.666667     6.052500
## 2          Cuyama  16.666667     3.624167
## 3   HeartofMojave   8.333333     5.637838
## 4    PanocheHills   8.333333     6.065833
## 5 SheepholeValley   0.000000     9.386916
## 6          Tecopa  20.833333     5.938333
## 7      TejonRanch  28.333333     2.535833
season1.frost <- season1 %>% group_by(Site)  %>% filter(year>2015) %>% summarize(count.consec(frost))
data.frame(season1.frost)
##              Site count.consec.frost.
## 1         Barstow                   5
## 2          Cuyama                   6
## 3   HeartofMojave                   4
## 4    PanocheHills                   2
## 5 SheepholeValley                   2
## 6          Tecopa                   4
## 7      TejonRanch                  10
season2.frost <- season2 %>% group_by(Site)  %>%  filter(year>2016)%>% summarize(count.consec(frost))
data.frame(season2.frost)
##              Site count.consec.frost.
## 1         Barstow                   5
## 2          Cuyama                   6
## 3   HeartofMojave                   3
## 4    PanocheHills                   3
## 5 SheepholeValley                   2
## 6          Tecopa                   5
## 7      TejonRanch                  10

Microenvironmental differences

We compared temperature and relative humidity between shrub and open microsites among all sites along the regional gradient


    Shapiro-Wilk normality test

data:  fit1$residuals
W = 0.87899, p-value < 2.2e-16
                     Df Sum Sq Mean Sq F value Pr(>F)    
Microsite             1   33.6   33.59  76.013 <2e-16 ***
gradient              6  234.0   39.00  88.248 <2e-16 ***
Microsite:gradient    6    8.8    1.46   3.306  0.003 ** 
Residuals          4290 1896.0    0.44                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness
Analysis of Deviance Table

Model: binomial, link: logit

Response: RH/100

Terms added sequentially (first to last)

                   Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
NULL                                4333     974.99             
Microsite           1     2.11      4332     972.88   0.1463    
gradient            6   528.28      4326     444.60   <2e-16 ***
Microsite:gradient  6     7.44      4320     437.17   0.2824    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Deviance Table

Model: Gamma, link: inverse

Response: swc

Terms added sequentially (first to last)

               Df Deviance Resid. Df Resid. Dev        F Pr(>F)    
NULL                             404    231.405                    
Site            6  189.683       398     41.722 381.2179 <2e-16 ***
Microsite       1    0.287       397     41.435   3.4579 0.0637 .  
Site:Microsite  6    0.696       391     40.739   1.3994 0.2136    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Deviance Table

Model: Gamma, link: inverse

Response: swc

Terms added sequentially (first to last)

               Df Deviance Resid. Df Resid. Dev        F Pr(>F)    
NULL                             419    183.193                    
Site            6  144.930       413     38.263 279.4537 <2e-16 ***
Microsite       1    0.080       412     38.183   0.9284 0.3358    
Site:Microsite  6    0.535       406     37.647   1.0322 0.4037    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

site level means

## site coordintes
site.gps <- read.csv("Data/ERGsites.csv")
# nutrient data for each site
nutrients <- read.csv("Data/ERG.soilnutrients.csv")
 ## long term climate data extracted from worldclim
clim.dat <- read.csv("Data/ERG.worldclim.csv")


## check correlation among long-term climate variables
cor(clim.dat[,3:8]) ## Temp.wettest.QR least correlated 
##                    Annual.Temp Temp.Seasonality Temp.wettest.QR
## Annual.Temp          1.0000000        0.9051791       0.7260617
## Temp.Seasonality     0.9051791        1.0000000       0.4650548
## Temp.wettest.QR      0.7260617        0.4650548       1.0000000
## Annual.precip       -0.8911069       -0.9650533      -0.5313460
## Precip.seasonality  -0.9715419       -0.9774646      -0.6165170
## Precipt.wettest.QR  -0.8839124       -0.9570957      -0.5523347
##                    Annual.precip Precip.seasonality Precipt.wettest.QR
## Annual.Temp           -0.8911069         -0.9715419         -0.8839124
## Temp.Seasonality      -0.9650533         -0.9774646         -0.9570957
## Temp.wettest.QR       -0.5313460         -0.6165170         -0.5523347
## Annual.precip          1.0000000          0.9613729          0.9980959
## Precip.seasonality     0.9613729          1.0000000          0.9543657
## Precipt.wettest.QR     0.9980959          0.9543657          1.0000000
## Obtain mean shrub traits for each site
shrubs <- read.csv("Data/ERG.shrub.csv")
shrubs <- subset(shrubs, Microsite=="shrub")
shrubs.mean <- shrubs %>% group_by(Gradient,Site) %>% summarise_all(funs(mean))
shrubs.mean <- data.frame(shrubs.mean)
shrubs.vars <- shrubs.mean[,c("Site","Gradient","volume","canopy","Dx","DxEph","Compaction")] 

## test correlation among shrub traits
cor(shrubs.vars[,3:7]) ## compaction x volume & DxEph x Dx high correlation
##                volume     canopy         Dx      DxEph Compaction
## volume      1.0000000  0.2243902  0.4003665  0.2631563 -0.7023062
## canopy      0.2243902  1.0000000 -0.6847957 -0.3538094 -0.1437590
## Dx          0.4003665 -0.6847957  1.0000000  0.7803527 -0.1355700
## DxEph       0.2631563 -0.3538094  0.7803527  1.0000000  0.3221944
## Compaction -0.7023062 -0.1437590 -0.1355700  0.3221944  1.0000000
vifstep(shrubs.vars[,3:7], th=10) ## Dx showing collinearity problems
## 1 variables from the 5 input variables have collinearity problem: 
##  
## Dx 
## 
## After excluding the collinear variables, the linear correlation coefficients ranges between: 
## min correlation ( Compaction ~ canopy ):  -0.143759 
## max correlation ( Compaction ~ volume ):  -0.7023062 
## 
## ---------- VIFs of the remained variables -------- 
##    Variables      VIF
## 1     volume 7.026144
## 2     canopy 1.925499
## 3      DxEph 4.316058
## 4 Compaction 6.400914
shrub.site <- shrubs.vars[,-c(1,2,5)]
rownames(shrub.site) <- shrubs.vars[,1]

## PCA of shrub characteristics
pca1 <- prcomp(log(shrub.site), scale=T)
plot(pca1)

biplot(pca1, scale = T)

summary(pca1) ## 77% variation explained
## Importance of components:
##                           PC1    PC2    PC3    PC4
## Standard deviation     1.3017 1.1778 0.9082 0.3059
## Proportion of Variance 0.4236 0.3468 0.2062 0.0234
## Cumulative Proportion  0.4236 0.7704 0.9766 1.0000
## PCA of site characteristics for season1
## Obtain mean weather variables for each site
season1.mean <- season1 %>% group_by(Gradient,Site) %>% summarise(temp.var=var(avg.temp,na.rm=T),Precip=sum(Precip),wind=mean(wind.speed, na.rm=T))
season1.mean  <- data.frame(season1.mean)
## extract key variables
## dropped RH, and chose min.temp because least correlation with others & cold stress

##combine nutrients and long-term averages
site.vars <- data.frame(season1.mean[,3:5],site.gps["elevation"]) ## drop Phosphorus and other climate variables because of correlations, 
row.names(site.vars) <- shrubs.vars[,1]
cor(site.vars)
##              temp.var     Precip        wind  elevation
## temp.var   1.00000000 -0.8679054  0.05816084 -0.2699588
## Precip    -0.86790537  1.0000000 -0.37279135  0.1222504
## wind       0.05816084 -0.3727914  1.00000000  0.1054497
## elevation -0.26995879  0.1222504  0.10544974  1.0000000
site.vars2016 <- site.vars

##  check for collinearity
vifstep(site.vars, th=10) ## remove potassium and temperature minimum
## No variable from the 4 input variables has collinearity problem. 
## 
## The linear correlation coefficients ranges between: 
## min correlation ( wind ~ temp.var ):  0.05816084 
## max correlation ( Precip ~ temp.var ):  -0.8679054 
## 
## ---------- VIFs of the remained variables -------- 
##   Variables      VIF
## 1  temp.var 6.640710
## 2    Precip 7.316232
## 3      wind 1.739558
## 4 elevation 1.142681
pca1 <- prcomp(log(abs(site.vars)), scale=T)
plot(pca1)

biplot(pca1)

## check contribution of loadings
pca1$rotation
##                  PC1        PC2         PC3        PC4
## temp.var   0.6001002 -0.1669018  0.50554765 -0.5970303
## Precip    -0.6703313 -0.1105034 -0.09936844 -0.7270288
## wind       0.2743589  0.7981791 -0.43410209 -0.3149488
## elevation -0.3395039  0.5681927  0.73898773  0.1256633
aload <- abs(pca1$rotation)
sweep(aload, 2, colSums(aload), "/")
##                 PC1       PC2        PC3        PC4
## temp.var  0.3184748 0.1015355 0.28433406 0.33832382
## Precip    0.3557466 0.0672253 0.05588758 0.41199111
## wind      0.1456030 0.4855763 0.24415109 0.17847449
## elevation 0.1801756 0.3456629 0.41562726 0.07121058
summary(pca1) ## 89% variation explained
## Importance of components:
##                          PC1    PC2    PC3     PC4
## Standard deviation     1.456 1.0458 0.8577 0.22479
## Proportion of Variance 0.530 0.2734 0.1839 0.01263
## Cumulative Proportion  0.530 0.8034 0.9874 1.00000
gradient1.season1 <- pca1$x[,1]
gradient2.season1 <- pca1$x[,2]

## season2
season2.mean <- season2 %>% group_by(Gradient,Site) %>% summarise(temp.var=var(avg.temp,na.rm=T),Precip=sum(Precip, na.rm=T),wind=mean(wind.speed, na.rm=T))
season2.mean  <- data.frame(season2.mean)
## extract key variables
## dropped RH, and chose min.temp because least correlation with others & cold stress

##combine nutrients and long-term averages
site.vars <- data.frame(season2.mean[,3:5],site.gps["elevation"]) ## drop Phosphorus and other climate variables because of correlations, 
row.names(site.vars) <-  shrubs.vars[,1]
cor(site.vars)
##             temp.var     Precip       wind  elevation
## temp.var   1.0000000 -0.7427028 -0.2262409 -0.3715471
## Precip    -0.7427028  1.0000000 -0.1551517  0.1813600
## wind      -0.2262409 -0.1551517  1.0000000  0.1052805
## elevation -0.3715471  0.1813600  0.1052805  1.0000000
site.vars2017 <- site.vars


##  check for collinearity
vifstep(site.vars, th=10) ## remove potassium and temperature minimum
## No variable from the 4 input variables has collinearity problem. 
## 
## The linear correlation coefficients ranges between: 
## min correlation ( elevation ~ wind ):  0.1052805 
## max correlation ( Precip ~ temp.var ):  -0.7427028 
## 
## ---------- VIFs of the remained variables -------- 
##   Variables      VIF
## 1  temp.var 3.438786
## 2    Precip 3.035337
## 3      wind 1.402010
## 4 elevation 1.192010
pca2 <- prcomp(log(abs(site.vars)), scale=T)
plot(pca2)

biplot(pca2)

summary(pca2) ## 85% variation explained
## Importance of components:
##                          PC1    PC2    PC3     PC4
## Standard deviation     1.434 1.0121 0.8809 0.37883
## Proportion of Variance 0.514 0.2561 0.1940 0.03588
## Cumulative Proportion  0.514 0.7701 0.9641 1.00000
## check contribution of loadings
pca2$rotation
##                  PC1         PC2         PC3       PC4
## temp.var   0.6557532 -0.09123357  0.20565681 0.7206729
## Precip    -0.6176537 -0.16891771 -0.40105332 0.6550778
## wind      -0.1019884  0.97279862  0.06826009 0.1964733
## elevation -0.4220070 -0.12963831  0.89005734 0.1135866
aload <- abs(pca2$rotation)
sweep(aload, 2, colSums(aload), "/")
##                  PC1        PC2       PC3        PC4
## temp.var  0.36483385 0.06695609 0.1314078 0.42749339
## Precip    0.34363685 0.12396828 0.2562596 0.38858328
## wind      0.05674213 0.71393442 0.0436159 0.11654531
## elevation 0.23478717 0.09514122 0.5687167 0.06737801
## define gradients
gradient1.season2 <- pca2$x[,1]
gradient2.season2 <- pca2$x[,2]

Compare gradient to phytometer

Differences from Nutrient content

## join aridity with nutrients
nutrients.climate <- merge(nutrients,subset(site.climate, Year==2016), by=c("Site"))


## Nitrogen difference between shrub and sites
m.nit <- aov(log(N) ~ Site * microsite, data=nutrients)
summary(m.nit)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## Site            6  63.04   10.51  13.241 2.86e-09 ***
## microsite       1  40.16   40.16  50.614 2.24e-09 ***
## Site:microsite  6   4.96    0.83   1.042    0.408    
## Residuals      56  44.43    0.79                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m.nit, "microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(N) ~ Site * microsite, data = nutrients)
## 
## $microsite
##                diff      lwr      upr p adj
## shrub-open 1.514873 1.088317 1.941429     0
## Potassium difference between shrub and sites
m.pot <- aov(log(K) ~ Site * microsite, data=nutrients)
summary(m.pot)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## Site            6 15.897   2.650  22.013 4.09e-13 ***
## microsite       1  4.445   4.445  36.934 1.14e-07 ***
## Site:microsite  6  1.605   0.268   2.223   0.0541 .  
## Residuals      56  6.740   0.120                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m.pot, "microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(K) ~ Site * microsite, data = nutrients)
## 
## $microsite
##                 diff       lwr       upr p adj
## shrub-open 0.5040081 0.3378755 0.6701406 1e-07
## Phosphorus difference between shrub and sites
m.pho <- aov(log(P) ~ Site * microsite, data=nutrients)
summary(m.pho)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## Site            6 23.163   3.861   14.33 8.10e-10 ***
## microsite       1 10.546  10.546   39.15 5.79e-08 ***
## Site:microsite  6  3.394   0.566    2.10   0.0676 .  
## Residuals      56 15.085   0.269                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m.pho, "microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(P) ~ Site * microsite, data = nutrients)
## 
## $microsite
##                 diff       lwr      upr p adj
## shrub-open 0.7762734 0.5277317 1.024815 1e-07
nutrient.mean <- nutrients.climate %>% group_by(Site, arid.gradient) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean)
 

ggplot(nutrient.mean) + geom_jitter(aes(x=arid.gradient, y=nitrogen), color="#E69F00", size=3) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=nitrogen), color="#E69F00") + geom_jitter(aes(x=arid.gradient, y=potassium/20), color="#56B4E9", size=3) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=potassium/20), color="#56B4E9") + geom_jitter(aes(x=arid.gradient, y=phosphorus), color="#009E73", size=3) +ylab("soil nutrient content")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("aridity gradient")+
annotate("segment", x = -4.5, xend = -4.35, y = 25, yend = 25,  colour = "#E69F00", size=2) + ## add custom legend
  annotate("text", x = -4.3, y = 25,  label = "Nitrogen", size=5, hjust=0)+
annotate("segment", x = -4.5, xend = -4.35, y = 23, yend = 23,  colour = "#56B4E9", size=2) + ## add custom legend
  annotate("text", x = -4.3, y = 23,  label = "Potassium", size=5, hjust=0)+
annotate("segment", x = -4.5, xend = -4.35, y = 21, yend = 21,  colour = "#009E73", size=2) + ## add custom legend
  annotate("text", x = -4.3, y = 21,  label = "Phosphorus", size=5, hjust=0)

m1 <- lm(nitrogen~poly(arid.gradient,2), data=nutrient.mean)
summary(m1)
## 
## Call:
## lm(formula = nitrogen ~ poly(arid.gradient, 2), data = nutrient.mean)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  1.1092  1.0570 -0.4565 -0.8511 -0.8789 -0.3757  0.3960 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               5.4326     0.3942  13.782 0.000161 ***
## poly(arid.gradient, 2)1 -12.6364     1.0429 -12.116 0.000266 ***
## poly(arid.gradient, 2)2   7.0832     1.0429   6.792 0.002454 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.043 on 4 degrees of freedom
## Multiple R-squared:  0.9797, Adjusted R-squared:  0.9695 
## F-statistic: 96.47 on 2 and 4 DF,  p-value: 0.0004126
m2 <- lm(phosphorus~arid.gradient, data=nutrient.mean)
summary(m2)
## 
## Call:
## lm(formula = phosphorus ~ arid.gradient, data = nutrient.mean)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  4.3035  1.1485 -2.0485  1.1930 -0.2791 -5.5020  1.1845 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     14.191      4.527   3.135   0.0258 *
## arid.gradient    1.780      1.298   1.371   0.2287  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.383 on 5 degrees of freedom
## Multiple R-squared:  0.2732, Adjusted R-squared:  0.1278 
## F-statistic: 1.879 on 1 and 5 DF,  p-value: 0.2287
m3 <- lm(potassium~poly(arid.gradient,2), data=nutrient.mean)
summary(m3)
## 
## Call:
## lm(formula = potassium ~ poly(arid.gradient, 2), data = nutrient.mean)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  -3.551 -65.927  32.611  57.492 -25.220  52.524 -47.931 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)   
## (Intercept)               186.56      22.69   8.221  0.00119 **
## poly(arid.gradient, 2)1   189.58      60.04   3.157  0.03427 * 
## poly(arid.gradient, 2)2   190.31      60.04   3.170  0.03387 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 60.04 on 4 degrees of freedom
## Multiple R-squared:  0.8334, Adjusted R-squared:  0.7502 
## F-statistic: 10.01 on 2 and 4 DF,  p-value: 0.02774

annual plant community repsonse

spp.data  <- read.csv("Data/ERG.communitydata.csv")
spp.data[is.na(spp.data)] <- 0

mean.spp <- spp.data %>% group_by(Year,Site, Microsite) %>%  summarize(abd=mean(Abundance), richness=mean(Richness), biomass=mean(Biomass))
mean.spp <- data.frame(mean.spp)

## collect aridity gradient data
site.climate <- rbind(s1.mean,s2.mean)
site.climate[,"Year"] <- c(rep("2016",7),rep("2017",7))

## combine climate data with community data
mean.spp <- merge(mean.spp, site.climate, by.y=c("Year","Site"))


## community models with GAMs
community <- spp.data
comm.clim <- merge(community,gams.data, by=c("Year","Site"))

m1 <- glmer.nb(Abundance~ Microsite * poly(gams,2) + (1|Year), data=comm.clim)
car::Anova(m1, type=c("II"))
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Abundance
##                             Chisq Df Pr(>Chisq)    
## Microsite                 19.0016  1  1.306e-05 ***
## poly(gams, 2)           1328.4369  2  < 2.2e-16 ***
## Microsite:poly(gams, 2)    3.1951  2     0.2024    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2 <- glmer.nb(Richness~ Microsite * poly(gams,2) + (1|Year), data=comm.clim)
## Warning in theta.ml(Y, mu, weights = object@resp$weights, limit = limit, :
## iteration limit reached
m2.poi <- glmer(Richness~ Microsite * poly(gams,2) + (1|Year), data=comm.clim, family="poisson")
car::Anova(m2.poi, type=c("III"))
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Richness
##                            Chisq Df Pr(>Chisq)    
## (Intercept)             135.1771  1  < 2.2e-16 ***
## Microsite                 7.8474  1   0.005089 ** 
## poly(gams, 2)            69.6946  2  7.345e-16 ***
## Microsite:poly(gams, 2)  32.9738  2  6.916e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 <- lmer(ihs(Biomass)~ Microsite * poly(gams,2) + (1|Year), data=comm.clim)
anova(m3)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                         Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Microsite                64.53  64.527     1 833.00  210.49 < 2.2e-16 ***
## poly(gams, 2)           511.61 255.805     2 833.44  834.46 < 2.2e-16 ***
## Microsite:poly(gams, 2)  32.08  16.041     2 833.00   52.33 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ## community responses season1
# mean.spp2016 <- subset(mean.spp, Year==2016)
# 
# ## richness
# plot(mean.spp2016[mean.spp2016$Microsite=="open","arid.gradient"], mean.spp2016[mean.spp2016$Microsite=="open","richness"], ylab="average site richness", xlab="aridity gradient")
# points(mean.spp2016[mean.spp2016$Microsite=="shrub","arid.gradient"],mean.spp2016[mean.spp2016$Microsite=="shrub","richness"], pch=19) 
# 
# m1 <- lm(richness~arid.gradient*Microsite, data=mean.spp2016)
# summary(m1)
# 
# ## abundance
# plot(mean.spp2016[mean.spp2016$Microsite=="open","arid.gradient"], mean.spp2016[mean.spp2016$Microsite=="open","abd"], ylab="average site abundance", xlab="aridity gradient")
# points(mean.spp2016[mean.spp2016$Microsite=="shrub","arid.gradient"],mean.spp2016[mean.spp2016$Microsite=="shrub","abd"], pch=19) 
# 
# m2 <- lm(abd~arid.gradient*Microsite, data=mean.spp2016)
# summary(m2)
# 
# ## biomass
# plot(mean.spp2016[mean.spp2016$Microsite=="open","arid.gradient"], mean.spp2016[mean.spp2016$Microsite=="open","biomass"], ylab="average site biomass", xlab="aridity gradient")
# points(mean.spp2016[mean.spp2016$Microsite=="shrub","arid.gradient"],mean.spp2016[mean.spp2016$Microsite=="shrub","biomass"], pch=19) 
# 
# m3 <- lm(biomass~arid.gradient*Microsite, data=mean.spp2016)
# summary(m3)
# 
# ##season2
# 
# ## community responses season2
# 
# mean.spp2017 <- subset(mean.spp, Year==2017)
# 
# 
# ## richness
# plot(mean.spp2017[mean.spp2017$Microsite=="open","arid.gradient"], mean.spp2017[mean.spp2017$Microsite=="open","richness"])
# points(mean.spp2017[mean.spp2017$Microsite=="shrub","arid.gradient"],mean.spp2017[mean.spp2017$Microsite=="shrub","richness"], pch=19) 
# 
# m1 <- lm(richness~arid.gradient*Microsite, data=mean.spp2017)
# summary(m1)
# 
# ## abundance
# plot(mean.spp2017[mean.spp2017$Microsite=="open","arid.gradient"], mean.spp2017[mean.spp2017$Microsite=="open","abd"])
# points(mean.spp2017[mean.spp2017$Microsite=="shrub","arid.gradient"],mean.spp2017[mean.spp2017$Microsite=="shrub","abd"], pch=19) 
# 
# m2 <- lm(abd~arid.gradient*Microsite, data=mean.spp2017)
# summary(m2)
# 
# ## biomass
# plot(mean.spp2017[mean.spp2017$Microsite=="open","arid.gradient"], mean.spp2017[mean.spp2017$Microsite=="open","biomass"])
# points(mean.spp2017[mean.spp2017$Microsite=="shrub","arid.gradient"],mean.spp2017[mean.spp2017$Microsite=="shrub","biomass"], pch=19) 
# 
# m3 <- lm(biomass~arid.gradient*Microsite, data=mean.spp2017)
# summary(m3)


# ## both seasons
# 
# ## richness
# plot(mean.spp[mean.spp$Microsite=="open","arid.gradient"], mean.spp[mean.spp$Microsite=="open","richness"],  ylab="average site richness", xlab="aridity gradient")
# points(mean.spp[mean.spp$Microsite=="shrub","arid.gradient"],mean.spp[mean.spp$Microsite=="shrub","richness"], pch=19) 
# 
# m1 <- lm(richness~arid.gradient*Microsite, data=mean.spp)
# summary(m1)
# 
# m1 <- lm(richness~poly(arid.gradient,2), data=mean.spp)
# summary(m1)
# 
# ggplot(mean.spp) + geom_jitter(aes(x=arid.gradient, y=richness, fill=Microsite, color=Microsite), size=3) + 
#  scale_color_manual(values=c("#E69F00","#56B4E9")) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=richness), color="black")+ ylab("Species richness") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 
# 
# ## abundance
# plot(mean.spp[mean.spp$Microsite=="open","arid.gradient"], mean.spp[mean.spp$Microsite=="open","abd"], ylab="average site abundance", xlab="aridity gradient")
# points(mean.spp[mean.spp$Microsite=="shrub","arid.gradient"],mean.spp[mean.spp$Microsite=="shrub","abd"], pch=19) 
# 
# m2 <- lm(abd~arid.gradient*Microsite, data=mean.spp)
# summary(m2)
# 
# m2 <- lm(ihs(abd)~arid.gradient, data=mean.spp)
# summary(m2)
# 
# ggplot(mean.spp) + geom_jitter(aes(x=arid.gradient, y=ihs(abd), fill=Microsite, color=Microsite), size=3)+ 
#  scale_color_manual(values=c("#E69F00","#56B4E9"))  +  stat_smooth(method="lm", formula= y~x,aes(x=arid.gradient, y=ihs(abd)), color="black")+ ylab("Annual abundance") + xlab("aridity gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 
# 
# 
# 
# ## biomass
# plot(mean.spp[mean.spp$Microsite=="open","arid.gradient"], ihs(mean.spp[mean.spp$Microsite=="open","biomass"]), ylab="average site biomass", xlab="aridity gradient")
# points(mean.spp[mean.spp$Microsite=="shrub","arid.gradient"],ihs(mean.spp[mean.spp$Microsite=="shrub","biomass"]), pch=19)
# 
# m3 <- lm(ihs(biomass)~arid.gradient*Microsite, data=mean.spp)
# summary(m3)
# 
# ggplot(mean.spp) + geom_jitter(aes(x=arid.gradient, y=ihs(biomass), fill=Microsite, color=Microsite), size=3)+ 
#  scale_color_manual(values=c("#E69F00","#56B4E9"))+  stat_smooth(method="lm", formula= y~x,aes(x=arid.gradient, y=ihs(biomass)), color="#56B4E9", data=subset(mean.spp, Microsite=="shrub"))  +  stat_smooth(method="lm", formula= y~x,aes(x=arid.gradient, y=ihs(biomass)), data=subset(mean.spp, Microsite=="open"), color="#E69F00")+ ylab("Annual biomass") + xlab("aridity gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 
# 

RDA community and environmental data

Phytometer analysis

Difference in effect size

# 
# library(bootES)
# 
# effect.cal <- function(x){ bootES(x, R=999, data.col="Biomass", group.col="Microsite", contrast=c(shrub=1,open=-1), effect.type=c("cohens.d"))}
# 
# effect.site <- by(subset(comm2016, Site != "Barstow"), comm2016$Site[comm2016$Site != "Barstow"], FUN=effect.cal)
# 
# site.summary <- rbind(summary(effect.site$Panoche),summary(effect.site$Cuyama),summary(effect.site$TejonRanch),data.frame(stat=0,ci.low=0,ci.high=0,bias=0,std.error=0),summary(effect.site$MojavePreserve),summary(effect.site$SheepholeValley),summary(effect.site$Tecopa)) 
# site.summary[5,] <- 0
# site.summary <- sapply(site.summary, as.numeric)
# 
# par(mar=c(4.5,4.5,.5,.5))
# plot(c(1:7),site.summary[,"stat"], xlim=c(0.5,7.5), ylim=c(-1,3), pch=19, cex=1.4, ylab="Hedge's G", xlab="", xaxt="n")
# arrows(c(1:7),as.numeric(site.summary[,"ci.high"]),c(1:7),site.summary[,"ci.low"], angle=90, code=3, length=0, lwd=2)
# abline(h=0, lty=2, lwd=2)
# axis(1, c(1:7), labels=site.names, cex.axis=1)
# 
# ## phytometer
# 
# effect.cal <- function(x){ bootES(x, R=999, data.col="phyto.biomass", group.col="Microsite", contrast=c(shrub=1,open=-1), effect.type=c("cohens.d"))}
# 
# census2016 <- subset(census, Census=="end" & Year ==2016)
# 
# effect.site <- by(subset(census2016, Site != "Barstow"), census2016$Site[census2016$Site != "Barstow"], FUN=effect.cal)
# 
# site.summary <- rbind(summary(effect.site$Panoche),summary(effect.site$Cuyama),summary(effect.site$TejonRanch),data.frame(stat=0,ci.low=0,ci.high=0,bias=0,std.error=0),summary(effect.site$MojavePreserve),summary(effect.site$SheepholeValley),summary(effect.site$Tecopa)) 
# site.summary <- sapply(site.summary, as.numeric)
# 
# par(mar=c(4.5,4.5,.5,.5))
# plot(c(1:7),site.summary[,"stat"], xlim=c(0.5,7.5), ylim=c(-1,3), pch=19, cex=1.4, ylab="Hedge's G", xlab="", xaxt="n")
# arrows(c(1:7),site.summary[,"ci.high"],c(1:7),site.summary[,"ci.low"], angle=90, code=3, length=0, lwd=2)
# abline(h=0, lty=2, lwd=2)
# axis(1, c(1:7), labels=site.names, cex.axis=1)

Rii between years

Native non-native comparisons

## load species list
spp.list <- read.csv("Data/ERG.specieslist.csv")

## convert data to long format
status <- gather(community, species, abundance, 13:53)
names(spp.list)[1] <- "species" ## rename species column for merge

## combine native-non-native status with data set
status <- merge(status, spp.list, by="species")

## summarize per plot native and non-natives, convert back to wide format
mean.status <- status %>% group_by(Year, Site, Microsite, Rep, status) %>% summarize(abundance=sum(abundance))
mean.status <- spread(mean.status, status, abundance)
status.data <- mean.status

## calculate site means for native and non-native
mean.status <- mean.status %>% group_by(Year, Site, Microsite) %>% summarize(native=mean(native), non.native=mean(non.native))

## merge with original community data
mean.spp <- merge(mean.spp, mean.status, by=c("Year","Site","Microsite"))

## native
m1 <- lm(ihs(native) ~ poly(arid.gradient,2), data=subset(mean.spp, Microsite=="shrub"))
summary(m1)
## 
## Call:
## lm(formula = ihs(native) ~ poly(arid.gradient, 2), data = subset(mean.spp, 
##     Microsite == "shrub"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.03332 -0.65782 -0.04555  0.60465  1.25478 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.6533     0.2185   7.565 1.11e-05 ***
## poly(arid.gradient, 2)1   0.3447     0.8177   0.422  0.68145    
## poly(arid.gradient, 2)2  -2.9468     0.8177  -3.604  0.00414 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8177 on 11 degrees of freedom
## Multiple R-squared:  0.5448, Adjusted R-squared:  0.4621 
## F-statistic: 6.583 on 2 and 11 DF,  p-value: 0.01318
m1 <- lm(ihs(native) ~ arid.gradient, data=subset(mean.spp, Microsite=="open"))
summary(m1)
## 
## Call:
## lm(formula = ihs(native) ~ arid.gradient, data = subset(mean.spp, 
##     Microsite == "open"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.89497 -0.73527  0.02181  0.50785  1.81716 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)     4.0097     1.0111   3.966  0.00187 **
## arid.gradient   0.6575     0.3294   1.996  0.06913 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.093 on 12 degrees of freedom
## Multiple R-squared:  0.2493, Adjusted R-squared:  0.1867 
## F-statistic: 3.984 on 1 and 12 DF,  p-value: 0.06913
m1 <- lm(ihs(native) ~ poly(arid.gradient,2)*Microsite, data=mean.spp)
summary(m1)
## 
## Call:
## lm(formula = ihs(native) ~ poly(arid.gradient, 2) * Microsite, 
##     data = mean.spp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6483 -0.6647 -0.1444  0.5383  2.0863 
## 
## Coefficients:
##                                        Estimate Std. Error t value
## (Intercept)                              2.0776     0.2616   7.943
## poly(arid.gradient, 2)1                  3.0858     1.3841   2.229
## poly(arid.gradient, 2)2                 -1.1145     1.3841  -0.805
## Micrositeshrub                          -0.4243     0.3699  -1.147
## poly(arid.gradient, 2)1:Micrositeshrub  -2.5983     1.9574  -1.327
## poly(arid.gradient, 2)2:Micrositeshrub  -3.0529     1.9574  -1.560
##                                        Pr(>|t|)    
## (Intercept)                            6.65e-08 ***
## poly(arid.gradient, 2)1                  0.0363 *  
## poly(arid.gradient, 2)2                  0.4293    
## Micrositeshrub                           0.2637    
## poly(arid.gradient, 2)1:Micrositeshrub   0.1980    
## poly(arid.gradient, 2)2:Micrositeshrub   0.1331    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9787 on 22 degrees of freedom
## Multiple R-squared:  0.4229, Adjusted R-squared:  0.2918 
## F-statistic: 3.225 on 5 and 22 DF,  p-value: 0.0247
ggplot(mean.spp) + geom_jitter(aes(x=arid.gradient, y=ihs(native), fill=Microsite, color=Microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))+  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=ihs(native)), color="black")+ ylab("Native abundance") + xlab("aridity gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

## non-native
m1 <- lm(ihs(non.native) ~ poly(arid.gradient,2), data=subset(mean.spp, Microsite=="shrub"))
summary(m1)
## 
## Call:
## lm(formula = ihs(non.native) ~ poly(arid.gradient, 2), data = subset(mean.spp, 
##     Microsite == "shrub"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.32477 -0.20348 -0.05303  0.53297  2.27555 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               3.0777     0.3340   9.214 1.66e-06 ***
## poly(arid.gradient, 2)1   6.5874     1.2498   5.271 0.000264 ***
## poly(arid.gradient, 2)2   0.7721     1.2498   0.618 0.549292    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.25 on 11 degrees of freedom
## Multiple R-squared:  0.7191, Adjusted R-squared:  0.668 
## F-statistic: 14.08 on 2 and 11 DF,  p-value: 0.0009268
m1 <- lm(ihs(non.native) ~ poly(arid.gradient,2), data=subset(mean.spp, Microsite=="open"))
summary(m1)
## 
## Call:
## lm(formula = ihs(non.native) ~ poly(arid.gradient, 2), data = subset(mean.spp, 
##     Microsite == "open"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.66301 -0.60211  0.01432  0.55858  1.91804 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               2.8844     0.2828  10.200 6.06e-07 ***
## poly(arid.gradient, 2)1   6.1497     1.0581   5.812 0.000117 ***
## poly(arid.gradient, 2)2   0.4615     1.0581   0.436 0.671141    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.058 on 11 degrees of freedom
## Multiple R-squared:  0.7554, Adjusted R-squared:  0.7109 
## F-statistic: 16.98 on 2 and 11 DF,  p-value: 0.0004331
m1 <- lm(ihs(non.native) ~ poly(arid.gradient,2)*Microsite, data=mean.spp)
summary(m1)
## 
## Call:
## lm(formula = ihs(non.native) ~ poly(arid.gradient, 2) * Microsite, 
##     data = mean.spp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.32477 -0.43001 -0.01113  0.58589  2.27555 
## 
## Coefficients:
##                                        Estimate Std. Error t value
## (Intercept)                              2.8844     0.3095   9.320
## poly(arid.gradient, 2)1                  8.6970     1.6376   5.311
## poly(arid.gradient, 2)2                  0.6527     1.6376   0.399
## Micrositeshrub                           0.1933     0.4377   0.442
## poly(arid.gradient, 2)1:Micrositeshrub   0.6190     2.3159   0.267
## poly(arid.gradient, 2)2:Micrositeshrub   0.4393     2.3159   0.190
##                                        Pr(>|t|)    
## (Intercept)                            4.28e-09 ***
## poly(arid.gradient, 2)1                2.49e-05 ***
## poly(arid.gradient, 2)2                   0.694    
## Micrositeshrub                            0.663    
## poly(arid.gradient, 2)1:Micrositeshrub    0.792    
## poly(arid.gradient, 2)2:Micrositeshrub    0.851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.158 on 22 degrees of freedom
## Multiple R-squared:  0.7361, Adjusted R-squared:  0.6761 
## F-statistic: 12.27 on 5 and 22 DF,  p-value: 9.2e-06
m1 <- lm(ihs(non.native) ~ arid.gradient*Microsite, data=mean.spp)
summary(m1)
## 
## Call:
## lm(formula = ihs(non.native) ~ arid.gradient * Microsite, data = mean.spp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5300 -0.4742  0.1206  0.7847  2.0747 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    8.3299     1.0394   8.014 3.06e-08 ***
## arid.gradient                  1.8530     0.3386   5.472 1.26e-05 ***
## Micrositeshrub                 0.5808     1.4699   0.395    0.696    
## arid.gradient:Micrositeshrub   0.1319     0.4789   0.275    0.785    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.124 on 24 degrees of freedom
## Multiple R-squared:  0.7289, Adjusted R-squared:  0.695 
## F-statistic: 21.51 on 3 and 24 DF,  p-value: 5.506e-07
ggplot(mean.spp) + geom_jitter(aes(x=arid.gradient, y=ihs(non.native), fill=Microsite, color=Microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))+  stat_smooth(method="lm", formula= y~x,aes(x=arid.gradient, y=ihs(non.native)), color="black")+ ylab("Non-native abundance") + xlab("aridity gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

## mixed models
m1 <- lmer(ihs(non.native) ~ arid.gradient*Microsite + (1|Year), data=mean.spp)
shapiro.test(resid(m1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m1)
## W = 0.96481, p-value = 0.4502
anova(m1)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                         Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## arid.gradient           81.053  81.053     1 23.956  73.255 9.498e-09 ***
## Microsite                0.197   0.197     1 23.000   0.178    0.6768    
## arid.gradient:Microsite  0.096   0.096     1 23.000   0.087    0.7712    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m1) ## R squared values 
##       R2m       R2c 
## 0.7152612 0.7820962
m2 <- lmer(ihs(native) ~ poly(arid.gradient,2)*Microsite + (1|Year), data=mean.spp)
shapiro.test(resid(m2))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2)
## W = 0.97665, p-value = 0.7641
anova(m2)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                                   Sum Sq Mean Sq NumDF DenDF F.value
## poly(arid.gradient, 2)           10.1669  5.0834     2    22  5.3071
## Microsite                         1.2604  1.2604     1    22  1.3159
## poly(arid.gradient, 2):Microsite  4.0179  2.0089     2    22  2.0973
##                                   Pr(>F)  
## poly(arid.gradient, 2)           0.01316 *
## Microsite                        0.26366  
## poly(arid.gradient, 2):Microsite 0.14666  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values 
##       R2m       R2c 
## 0.3739082 0.3739082

Phylogenetics

library(picante)
library(ape)
library(brranching)
library(taxize)

## load species list to create tree
spp.list <- read.csv("Data/ERG.specieslist.csv")

## create a combined species name column and drop erinoginum because not to species level
taxa <- paste(spp.list$Genus, spp.list$Species.name)
taxa <- taxa[taxa!="Erinoginum spp"]

## create tree of phylogeny
tree <- phylomatic(taxa=taxa, get = 'POST')

## calculate branch lengths using  Grafen's method
## Grafen, A. (1989) The phylogenetic regression. Philosophical Transactions of the Royal society of London. Series B. Biological Sciences, 326, 119–157.
tree2 <- compute.brlen(tree, method="Grafen")

## view tree and outputs to verify branches
plot(tree2)

tree2
## 
## Phylogenetic tree with 39 tips and 29 internal nodes.
## 
## Tip labels:
##  cryptantha_barbigera, cryptantha_intermedia, phacelia_crenulata, phacelia_tanacetifolia, pectocarya_spp, amsinckia_grandiflora, ...
## Node labels:
##  poales_to_asterales, eudicots, , , asterids, ericales_to_asterales, ...
## 
## Rooted; includes branch lengths.
## format community to only have site and microsite
comm <- community %>% group_by(Year, Site, Microsite) %>% summarise_if(is.numeric, funs(sum(., na.rm=T)))
comm <- comm[,c(1:3,13:53)] ## extract species abundances only
comm <- data.frame(comm) ## convert to dataframe
comm[,"site.micro"] <- paste(comm$Site,comm$Microsite,as.character(comm$Year)) ## create site by microsite column
rownames(comm) <- comm[,length(comm)] ## replace row names with site by microsite
comm[,c("Year","Site","Microsite","Erinoginum.spp","Boraginaceae.sp.","site.micro")] <- NULL ## drop all columns but ones for analysis

## match species name format
spp.list[,"spp.name"] <- paste(spp.list$Genus, spp.list$Species.name)
colnames(comm) <- spp.list$spp.name[match(colnames(comm), spp.list$Species.shorthand)]

## match species format
names(comm) <- gsub(" ", "_", names(comm), fixed=T) ## underscores instead of spaces
names(comm) <- gsub("__", "_", names(comm), fixed=T) ## replace double underscores with one
names(comm) <- tolower(names(comm)) ## lower case names


## mean phylogenetic distance
mpd.data <- mpd(comm, cophenetic(tree2), abundance.weighted=T)

## format dataframe
mpd.data <- data.frame(Year=c(rep(2016,14),rep(2017,14)),Microsite=c("open","shrub"),site=rownames(comm), mpd=mpd.data)
mpd.data[,"Site"] <- gsub(" open", "", mpd.data[,"site"]) 
mpd.data[,"Site"] <- gsub(" shrub", "", mpd.data[,"Site"]) 
mpd.data[,"Site"] <- gsub(" 2016", "", mpd.data[,"Site"]) 
mpd.data[,"Site"] <- gsub(" 2017", "", mpd.data[,"Site"]) 
mpd.data[,"site"] <- NULL
mpd.data[is.na(mpd.data)] <- 0 ## barstow no plants
 

## getspecies richness for sites
spp.rich <- community %>% group_by(Year, Site, Microsite) %>% summarize(richness=mean(Richness))
spp.rich <- data.frame(spp.rich)

## join data
mpd.rich <- merge(mpd.data, spp.rich, by=c("Site","Microsite","Year"))

mpd.arid <- merge(mpd.rich, mean.spp, by=c("Site","Microsite","Year"))


## phylogenetic distance
mpd.arid <- mpd.arid[c(2,4:nrow(mpd.arid)),] ## drop barstow

m1 <- lm(mpd~arid.gradient*Microsite, data=mpd.arid)
summary(m1)
## 
## Call:
## lm(formula = mpd ~ arid.gradient * Microsite, data = mpd.arid)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.51050 -0.19590  0.06299  0.14972  0.71607 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                   0.97618    0.30321   3.219  0.00395 **
## arid.gradient                 0.04716    0.10394   0.454  0.65451   
## Micrositeshrub               -1.07891    0.42881  -2.516  0.01966 * 
## arid.gradient:Micrositeshrub -0.30441    0.14700  -2.071  0.05031 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2938 on 22 degrees of freedom
## Multiple R-squared:  0.3146, Adjusted R-squared:  0.2211 
## F-statistic: 3.365 on 3 and 22 DF,  p-value: 0.03692
ggplot(mpd.arid) + geom_jitter(aes(x=arid.gradient, y=mpd, fill=Microsite, color=Microsite), size=3) + 
 scale_color_manual(values=c("#E69F00","#56B4E9")) +  stat_smooth(method="lm", formula= y~x,aes(x=arid.gradient, y=mpd), data=subset(mpd.arid, Microsite=="shrub"), color="#56B4E9")+  stat_smooth(method="lm", formula= y~x,aes(x=arid.gradient, y=mpd), data=subset(mpd.arid, Microsite=="open"), color="#E69F00")+ ylab("Phylogenetic Community Dissimilarity") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

m1 <- lm(mpd~arid.gradient, data=subset(mpd.arid, Microsite=="shrub"))
summary(m1)
## 
## Call:
## lm(formula = mpd ~ arid.gradient, data = subset(mpd.arid, Microsite == 
##     "shrub"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.51050 -0.23112  0.03664  0.14204  0.71607 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    -0.1027     0.3563  -0.288    0.778  
## arid.gradient  -0.2573     0.1222  -2.106    0.059 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3452 on 11 degrees of freedom
## Multiple R-squared:  0.2873, Adjusted R-squared:  0.2225 
## F-statistic: 4.435 on 1 and 11 DF,  p-value: 0.05899
m2 <- lm(mpd~arid.gradient, data=subset(mpd.arid, Microsite=="open"))
summary(m2)
## 
## Call:
## lm(formula = mpd ~ arid.gradient, data = subset(mpd.arid, Microsite == 
##     "open"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38078 -0.09025  0.08952  0.15668  0.33783 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    0.97618    0.23852   4.093  0.00178 **
## arid.gradient  0.04716    0.08177   0.577  0.57575   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2311 on 11 degrees of freedom
## Multiple R-squared:  0.02935,    Adjusted R-squared:  -0.05889 
## F-statistic: 0.3326 on 1 and 11 DF,  p-value: 0.5758

Changes in Mechanisms along gradient

### Nutrients 
nutrient.mean <- nutrients.climate %>% group_by(Site, arid.gradient, microsite) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean)

## nitrogen
m1 <- lm(nitrogen ~ arid.gradient, data=subset(nutrient.mean, microsite=="open"))
summary(m1)
## 
## Call:
## lm(formula = nitrogen ~ arid.gradient, data = subset(nutrient.mean, 
##     microsite == "open"))
## 
## Residuals:
##        1        3        5        7        9       11       13 
##  1.62742  0.10479 -1.35216  0.97441 -0.05162 -0.76608 -0.53677 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    -2.7141     1.5041  -1.804   0.1310  
## arid.gradient  -1.2879     0.4313  -2.986   0.0306 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.124 on 5 degrees of freedom
## Multiple R-squared:  0.6407, Adjusted R-squared:  0.5689 
## F-statistic: 8.917 on 1 and 5 DF,  p-value: 0.03058
anova(m1)
## Analysis of Variance Table
## 
## Response: nitrogen
##               Df Sum Sq Mean Sq F value  Pr(>F)  
## arid.gradient  1 11.263  11.263  8.9174 0.03058 *
## Residuals      5  6.315   1.263                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(nutrient.mean) + geom_jitter(aes(x=arid.gradient, y=nitrogen, fill=microsite, color=microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9")) +  
  stat_smooth(method="lm", formula= y~x,aes(x=arid.gradient, y=nitrogen), color="#56B4E9", data=subset(nutrient.mean, microsite=="shrub")) +  
  stat_smooth(method="lm", formula= y~x,aes(x=arid.gradient, y=nitrogen), color="#E69F00", data=subset(nutrient.mean, microsite=="open"))+ ylab("Nitrogen") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

## phosphorus
m2 <- lm(phosphorus ~ arid.gradient * microsite, data=nutrient.mean)
summary(m2)
## 
## Call:
## lm(formula = phosphorus ~ arid.gradient * microsite, data = nutrient.mean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1684 -2.3472  0.6103  1.9735  7.0340 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                    8.8933     5.3802   1.653    0.129
## arid.gradient                  0.7555     1.5428   0.490    0.635
## micrositeshrub                10.5952     7.6088   1.392    0.194
## arid.gradient:micrositeshrub   2.0483     2.1819   0.939    0.370
## 
## Residual standard error: 4.02 on 10 degrees of freedom
## Multiple R-squared:  0.3967, Adjusted R-squared:  0.2158 
## F-statistic: 2.192 on 3 and 10 DF,  p-value: 0.152
anova(m2)
## Analysis of Variance Table
## 
## Response: phosphorus
##                         Df  Sum Sq Mean Sq F value Pr(>F)
## arid.gradient            1  43.011  43.011  2.6614 0.1339
## microsite                1  49.031  49.031  3.0339 0.1122
## arid.gradient:microsite  1  14.243  14.243  0.8813 0.3700
## Residuals               10 161.610  16.161
ggplot(nutrient.mean) + geom_jitter(aes(x=arid.gradient, y=phosphorus, fill=microsite, color=microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))  ## not significant

m3 <- lm(potassium ~ arid.gradient * microsite, data=nutrient.mean)
summary(m3)
## 
## Call:
## lm(formula = potassium ~ arid.gradient * microsite, data = nutrient.mean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -160.67  -64.64  -14.22   64.63  206.92 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                    290.78     145.57   1.998   0.0737 .
## arid.gradient                   44.64      41.74   1.069   0.3101  
## micrositeshrub                 278.33     205.87   1.352   0.2062  
## arid.gradient:micrositeshrub    56.24      59.03   0.953   0.3633  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 108.8 on 10 degrees of freedom
## Multiple R-squared:  0.4843, Adjusted R-squared:  0.3295 
## F-statistic:  3.13 on 3 and 10 DF,  p-value: 0.07444
anova(m3)
## Analysis of Variance Table
## 
## Response: potassium
##                         Df Sum Sq Mean Sq F value Pr(>F)  
## arid.gradient            1  71879   71879  6.0754 0.0334 *
## microsite                1  28476   28476  2.4069 0.1519  
## arid.gradient:microsite  1  10736   10736  0.9074 0.3633  
## Residuals               10 118313   11831                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## potassium
ggplot(nutrient.mean) + geom_jitter(aes(x=arid.gradient, y=potassium, fill=microsite, color=microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9")) +  
  stat_smooth(method="lm", formula= y~x,aes(x=arid.gradient, y=potassium), color="black") 

#### Shrubs increase nitrogen availability at aridity extremes #### Shrubs do not effect potassium or phosphorus availability

se <- function(x) sd(na.omit(x))/sqrt(length(na.omit(x)))


smc <- census %>% group_by(Year, Site, Microsite, Census) %>% summarize(swc.avg=mean(swc), swc.se=se(swc))

smc.arid <- merge(smc, site.climate, by=c("Site","Year"))


## early swc
smc.early <- subset(smc.arid, Census=="emergence")

ggplot(smc.early)  + geom_jitter(aes(x=arid.gradient, y=swc.avg, fill=Microsite, color=Microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))+  
  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=swc.avg), color="black") + ylab("Percent soil Moisture Content (emergence)") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

m1 <- lm(swc.avg~ poly(arid.gradient,2) * Microsite, data=smc.early)
summary(m1)
## 
## Call:
## lm(formula = swc.avg ~ poly(arid.gradient, 2) * Microsite, data = smc.early)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8989 -2.2340  0.1864  2.2093  8.7293 
## 
## Coefficients:
##                                        Estimate Std. Error t value
## (Intercept)                             11.1650     1.3550   8.240
## poly(arid.gradient, 2)1                 32.2475     7.1700   4.498
## poly(arid.gradient, 2)2                 20.8675     7.1700   2.910
## Micrositeshrub                          -0.4679     1.9163  -0.244
## poly(arid.gradient, 2)1:Micrositeshrub  -3.3904    10.1399  -0.334
## poly(arid.gradient, 2)2:Micrositeshrub  -3.1002    10.1399  -0.306
##                                        Pr(>|t|)    
## (Intercept)                             3.6e-08 ***
## poly(arid.gradient, 2)1                0.000179 ***
## poly(arid.gradient, 2)2                0.008109 ** 
## Micrositeshrub                         0.809376    
## poly(arid.gradient, 2)1:Micrositeshrub 0.741276    
## poly(arid.gradient, 2)2:Micrositeshrub 0.762672    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.07 on 22 degrees of freedom
## Multiple R-squared:  0.699,  Adjusted R-squared:  0.6306 
## F-statistic: 10.22 on 5 and 22 DF,  p-value: 3.647e-05
## end swc
smc.end <- subset(smc.arid, Census=="end")

ggplot(smc.end)  + geom_jitter(aes(x=arid.gradient, y=swc.avg, fill=Microsite, color=Microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))+  
  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=swc.avg), color="black") + ylab("Percent soil Moisture Content (end)") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

m1 <- lm(swc.avg~ poly(arid.gradient,2) * Microsite, data=smc.end)
summary(m1)
## 
## Call:
## lm(formula = swc.avg ~ poly(arid.gradient, 2) * Microsite, data = smc.end)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7611 -2.4087  0.0009  2.0666  9.3760 
## 
## Coefficients:
##                                        Estimate Std. Error t value
## (Intercept)                               5.093      1.335   3.815
## poly(arid.gradient, 2)1                  15.502      7.065   2.194
## poly(arid.gradient, 2)2                  17.571      7.065   2.487
## Micrositeshrub                           -0.576      1.888  -0.305
## poly(arid.gradient, 2)1:Micrositeshrub   -1.945      9.991  -0.195
## poly(arid.gradient, 2)2:Micrositeshrub   -3.453      9.991  -0.346
##                                        Pr(>|t|)    
## (Intercept)                            0.000946 ***
## poly(arid.gradient, 2)1                0.039073 *  
## poly(arid.gradient, 2)2                0.020948 *  
## Micrositeshrub                         0.763211    
## poly(arid.gradient, 2)1:Micrositeshrub 0.847426    
## poly(arid.gradient, 2)2:Micrositeshrub 0.732914    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.996 on 22 degrees of freedom
## Multiple R-squared:  0.4604, Adjusted R-squared:  0.3377 
## F-statistic: 3.754 on 5 and 22 DF,  p-value: 0.01312

Shrubs do not alter soil moisture availability

hobo.means <- HOBO.data %>% group_by(season, Microsite, Site) %>% summarize(temp.avg=mean(Temp),temp.var=var(Temp),temp.se=se(Temp),rh.avg=mean(RH, na.rm=T),rh.var=var(RH, na.rm=T),rh.se=se(RH),frost=(sum(na.omit(Temp)<0)/length(na.omit(Temp))*100),heat=(sum(na.omit(Temp)>30)/length(na.omit(Temp))*100))
hobo.means[,"Year"] <- c(rep(2016,14),rep(2017,14))


hobo.arid <- merge(hobo.means, site.climate, by=c("Site","Year"))

## temperature
ggplot(hobo.arid)  + geom_jitter(aes(x=arid.gradient, y=temp.avg, fill=Microsite, color=Microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))+  
  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=temp.avg), color="black") + ylab("Average Temperature (C°)") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

m1 <- lm(temp.avg~arid.gradient * Microsite, data=hobo.arid)
summary(m1)
## 
## Call:
## lm(formula = temp.avg ~ arid.gradient * Microsite, data = hobo.arid)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2362 -0.7204 -0.3884  0.4197  3.2981 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   6.21950    1.38277   4.498 0.000149 ***
## arid.gradient                -1.75909    0.45045  -3.905 0.000669 ***
## Micrositeshrub               -1.42763    1.95554  -0.730 0.472428    
## arid.gradient:Micrositeshrub  0.00911    0.63704   0.014 0.988709    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.495 on 24 degrees of freedom
## Multiple R-squared:  0.6064, Adjusted R-squared:  0.5571 
## F-statistic: 12.32 on 3 and 24 DF,  p-value: 4.451e-05
anova(m1)
## Analysis of Variance Table
## 
## Response: temp.avg
##                         Df Sum Sq Mean Sq F value    Pr(>F)    
## arid.gradient            1 67.815  67.815 30.3425 1.154e-05 ***
## Microsite                1 14.807  14.807  6.6252   0.01666 *  
## arid.gradient:Microsite  1  0.000   0.000  0.0002   0.98871    
## Residuals               24 53.639   2.235                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.s <- lm(temp.avg~arid.gradient, data=subset(hobo.arid, Microsite=="shrub"))
summary(m1.s)
## 
## Call:
## lm(formula = temp.avg ~ arid.gradient, data = subset(hobo.arid, 
##     Microsite == "shrub"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5836 -0.9458 -0.4008  0.4122  3.2981 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)     4.7919     1.3727   3.491  0.00446 **
## arid.gradient  -1.7500     0.4472  -3.913  0.00206 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.484 on 12 degrees of freedom
## Multiple R-squared:  0.5607, Adjusted R-squared:  0.5241 
## F-statistic: 15.31 on 1 and 12 DF,  p-value: 0.00206
## RH
ggplot(hobo.arid)  + geom_jitter(aes(x=arid.gradient, y=rh.avg, fill=Microsite, color=Microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))+  
  stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=rh.avg), color="black") + ylab("Relative Humidity (%)") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

m2 <- lm(log(rh.avg)~arid.gradient * Microsite, data=hobo.arid)
summary(m2)
## 
## Call:
## lm(formula = log(rh.avg) ~ arid.gradient * Microsite, data = hobo.arid)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.313369 -0.094086  0.008993  0.111107  0.252701 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  4.703640   0.165465  28.427  < 2e-16 ***
## arid.gradient                0.251051   0.053902   4.658 9.92e-05 ***
## Micrositeshrub               0.049072   0.234003   0.210    0.836    
## arid.gradient:Micrositeshrub 0.004987   0.076229   0.065    0.948    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1789 on 24 degrees of freedom
## Multiple R-squared:  0.6497, Adjusted R-squared:  0.6059 
## F-statistic: 14.84 on 3 and 24 DF,  p-value: 1.131e-05
anova(m2)
## Analysis of Variance Table
## 
## Response: log(rh.avg)
##                         Df  Sum Sq Mean Sq F value    Pr(>F)    
## arid.gradient            1 1.41614 1.41614 44.2510 7.014e-07 ***
## Microsite                1 0.00829 0.00829  0.2591    0.6154    
## arid.gradient:Microsite  1 0.00014 0.00014  0.0043    0.9484    
## Residuals               24 0.76806 0.03200                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## frost
hobo.arid2 <- subset(hobo.arid, Site!="Cuyama") ## drop cuyama because abnormal frost periods

ggplot(hobo.arid2)  + geom_jitter(aes(x=arid.gradient, y=frost, fill=Microsite, color=Microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))+  
  stat_smooth(method="lm", formula= y~x,aes(x=arid.gradient, y=frost), color="black") + ylab("Number of frost days") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

m3 <- lm(frost~arid.gradient * Microsite, data=hobo.arid2)
summary(m3)
## 
## Call:
## lm(formula = frost ~ arid.gradient * Microsite, data = hobo.arid2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2695 -1.3394 -0.3817  1.9954  3.3744 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   10.7896     2.2613   4.771 0.000116 ***
## arid.gradient                  1.4999     0.7049   2.128 0.045970 *  
## Micrositeshrub                -2.3006     3.1980  -0.719 0.480218    
## arid.gradient:Micrositeshrub  -0.3756     0.9969  -0.377 0.710319    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.117 on 20 degrees of freedom
## Multiple R-squared:  0.3059, Adjusted R-squared:  0.2018 
## F-statistic: 2.938 on 3 and 20 DF,  p-value: 0.05818
anova(m3)
## Analysis of Variance Table
## 
## Response: frost
##                         Df Sum Sq Mean Sq F value  Pr(>F)  
## arid.gradient            1 31.066 31.0655  6.9304 0.01596 *
## Microsite                1  7.806  7.8058  1.7414 0.20187  
## arid.gradient:Microsite  1  0.636  0.6363  0.1419 0.71032  
## Residuals               20 89.649  4.4825                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## temperature variability

ggplot(hobo.arid)  + geom_jitter(aes(x=arid.gradient, y=temp.var.x, fill=Microsite, color=Microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))+  
  stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=temp.var.x), color="#E69F00", data=subset(hobo.arid, Microsite=="open")) + ylab("Number of frost days")+  
  stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=temp.var.x), color="#56B4E9", data=subset(hobo.arid, Microsite=="shrub")) + ylab("temperature variation") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

m4 <- lm(log(temp.var.x)~arid.gradient * Microsite, data=hobo.arid)
summary(m4)
## 
## Call:
## lm(formula = log(temp.var.x) ~ arid.gradient * Microsite, data = hobo.arid)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.75540 -0.15564  0.02433  0.17182  0.46013 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   3.83222    0.26061  14.705 1.67e-13 ***
## arid.gradient                -0.19133    0.08490  -2.254   0.0336 *  
## Micrositeshrub               -0.46416    0.36856  -1.259   0.2200    
## arid.gradient:Micrositeshrub -0.02535    0.12006  -0.211   0.8346    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2818 on 24 degrees of freedom
## Multiple R-squared:   0.51,  Adjusted R-squared:  0.4488 
## F-statistic: 8.327 on 3 and 24 DF,  p-value: 0.0005707
anova(m4)
## Analysis of Variance Table
## 
## Response: log(temp.var.x)
##                         Df  Sum Sq Mean Sq F value   Pr(>F)   
## arid.gradient            1 0.91677 0.91677 11.5480 0.002368 **
## Microsite                1 1.06294 1.06294 13.3892 0.001241 **
## arid.gradient:Microsite  1 0.00354 0.00354  0.0446 0.834588   
## Residuals               24 1.90531 0.07939                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##  extreme heat
## reverse inverse hyperbolic sine
hs <- function(x) {
    y <- 0.5*exp(-x)*(exp(2*x)-1)
    return(y)
}

ggplot(hobo.arid)  + geom_jitter(aes(x=arid.gradient, y=heat, fill=Microsite, color=Microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))+  
  stat_smooth(method="lm", formula= y~hs(x),aes(x=arid.gradient, y=heat), color="#E69F00", data=subset(hobo.arid, Microsite=="open")) +  
  stat_smooth(method="lm", formula= y~hs(x),aes(x=arid.gradient, y=heat), color="#56B4E9", data=subset(hobo.arid, Microsite=="shrub")) + ylab("extreme heat days variation") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

m5 <- lm(ihs(heat)~arid.gradient * Microsite, data=hobo.arid)
summary(m5)
## 
## Call:
## lm(formula = ihs(heat) ~ arid.gradient * Microsite, data = hobo.arid)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7532 -0.3832 -0.1102  0.2072  1.2432 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -0.2199     0.4899  -0.449 0.657489    
## arid.gradient                 -0.7396     0.1596  -4.635 0.000105 ***
## Micrositeshrub                -0.1480     0.6928  -0.214 0.832615    
## arid.gradient:Micrositeshrub   0.1787     0.2257   0.792 0.436163    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5296 on 24 degrees of freedom
## Multiple R-squared:  0.6529, Adjusted R-squared:  0.6095 
## F-statistic: 15.05 on 3 and 24 DF,  p-value: 1.015e-05
anova(m5)
## Analysis of Variance Table
## 
## Response: ihs(heat)
##                         Df Sum Sq Mean Sq F value    Pr(>F)    
## arid.gradient            1 9.3144  9.3144 33.2088 6.124e-06 ***
## Microsite                1 3.1726  3.1726 11.3114   0.00258 ** 
## arid.gradient:Microsite  1 0.1759  0.1759  0.6271   0.43616    
## Residuals               24 6.7315  0.2805                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

shrubs reduce temperature variability

shrubz <- read.csv("Data/ERG.shrub.csv")

shrub.mean <- shrubz %>% group_by(Site, Microsite) %>% summarize(compact=mean(Compaction))

shrub.clim <- merge(shrub.mean, site.climate, by=c("Site"))
#shrub.clim <- subset(shrub.clim, Year==2016)

ggplot(shrub.clim)  + geom_jitter(aes(x=arid.gradient, y=compact, fill=Microsite, color=Microsite), size=3)+ 
 scale_color_manual(values=c("#E69F00","#56B4E9"))+   stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=compact), color="#E69F00", data=subset(shrub.clim, Microsite=="open"))+
  stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=compact), color="#56B4E9", data=subset(shrub.clim, Microsite=="shrub")) + ylab("soil compaction") + xlab("Aridity Gradient")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) 

m1 <- lm(log(compact) ~ arid.gradient* Microsite, data=shrub.clim)
summary(m1)
## 
## Call:
## lm(formula = log(compact) ~ arid.gradient * Microsite, data = shrub.clim)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.36122 -0.13659 -0.02823  0.05479  0.61368 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                   0.26645    0.23301   1.144  0.26410   
## arid.gradient                -0.07328    0.07590  -0.965  0.34394   
## Micrositeshrub               -1.19983    0.32952  -3.641  0.00130 **
## arid.gradient:Micrositeshrub -0.32933    0.10734  -3.068  0.00528 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2519 on 24 degrees of freedom
## Multiple R-squared:  0.5932, Adjusted R-squared:  0.5424 
## F-statistic: 11.67 on 3 and 24 DF,  p-value: 6.53e-05
anova(m1)
## Analysis of Variance Table
## 
## Response: log(compact)
##                         Df  Sum Sq Mean Sq F value    Pr(>F)    
## arid.gradient            1 1.24730 1.24730 19.6548 0.0001755 ***
## Microsite                1 0.37672 0.37672  5.9362 0.0226248 *  
## arid.gradient:Microsite  1 0.59733 0.59733  9.4127 0.0052769 ** 
## Residuals               24 1.52305 0.06346                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

shrubs reduce soil compaction and increase substrate penetration

Phytometer models

library(lsmeans)
## The 'lsmeans' package is being deprecated.
## Users are encouraged to switch to 'emmeans'.
## See help('transition') for more information, including how
## to convert 'lsmeans' objects and scripts to work with 'emmeans'.
## 
## Attaching package: 'lsmeans'
## The following object is masked from 'package:lmerTest':
## 
##     lsmeans
census[is.na(census)] <- 0
census[,"phyto.abd"] <- rowSums(census[,c("Phacelia","Plantago","Salvia")])

census.end <- subset(census, Census=="end")
census.end[,"Year"] <- as.factor(census.end$Year)
census.arid <- merge(census.end, site.climate, by=c("Year","Site"))


## all species end of season
m1 <- glm.nb(phyto.abd~ Microsite * Nutrient * arid.gradient, data=census.arid)
anova(m2, test="Chisq") ## Site and microsite*nutrient significant
## Analysis of Variance Table
## 
## Response: log(rh.avg)
##                         Df  Sum Sq Mean Sq F value    Pr(>F)    
## arid.gradient            1 1.41614 1.41614 44.2510 7.014e-07 ***
## Microsite                1 0.00829 0.00829  0.2591    0.6154    
## arid.gradient:Microsite  1 0.00014 0.00014  0.0043    0.9484    
## Residuals               24 0.76806 0.03200                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## run full model for phacelia biomass
m1 <- aov(log(Phacelia.biomass)~ Microsite * Nutrient * arid.gradient,  data=subset(census.arid, Phacelia.biomass>0))
summary(m1) ## site and microsite
##                                   Df Sum Sq Mean Sq F value   Pr(>F)    
## Microsite                          1   19.0   19.00   9.126  0.00281 ** 
## Nutrient                           1    8.8    8.76   4.207  0.04141 *  
## arid.gradient                      1  122.7  122.66  58.927 4.95e-13 ***
## Microsite:Nutrient                 1    0.2    0.16   0.078  0.78078    
## Microsite:arid.gradient            1    1.4    1.41   0.678  0.41121    
## Nutrient:arid.gradient             1    0.0    0.00   0.000  0.99128    
## Microsite:Nutrient:arid.gradient   1    2.3    2.27   1.088  0.29795    
## Residuals                        225  468.3    2.08                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m1, "Microsite")
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Microsite, arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Nutrient, arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Microsite, Nutrient, arid.gradient
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Phacelia.biomass) ~ Microsite * Nutrient * arid.gradient, data = subset(census.arid, Phacelia.biomass > 0))
## 
## $Microsite
##                 diff       lwr       upr     p adj
## shrub-open 0.5995862 0.2084783 0.9906941 0.0028111
TukeyHSD(m1, "Nutrient")
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Microsite, arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Nutrient, arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Microsite, Nutrient, arid.gradient
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Phacelia.biomass) ~ Microsite * Nutrient * arid.gradient, data = subset(census.arid, Phacelia.biomass > 0))
## 
## $Nutrient
##                     diff        lwr         upr     p adj
## Ambient-Added -0.3865601 -0.7590709 -0.01404918 0.0420304
## run full model for plantago biomass
m2 <- aov(log(Plantago.biomass)~ Microsite * Nutrient * arid.gradient,  data=subset(census.arid, Plantago.biomass>0))
summary(m2) ## site and year + sitexmicrosite
##                                   Df Sum Sq Mean Sq F value   Pr(>F)    
## Microsite                          1   9.59    9.59   5.718 0.018067 *  
## Nutrient                           1  26.14   26.14  15.582 0.000122 ***
## arid.gradient                      1 136.39  136.39  81.293 1.01e-15 ***
## Microsite:Nutrient                 1   0.61    0.61   0.366 0.546337    
## Microsite:arid.gradient            1   1.57    1.57   0.937 0.334540    
## Nutrient:arid.gradient             1   3.96    3.96   2.362 0.126491    
## Microsite:Nutrient:arid.gradient   1   3.21    3.21   1.912 0.168873    
## Residuals                        146 244.96    1.68                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m2, "Microsite")
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Microsite, arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Nutrient, arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Microsite, Nutrient, arid.gradient
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Plantago.biomass) ~ Microsite * Nutrient * arid.gradient, data = subset(census.arid, Plantago.biomass > 0))
## 
## $Microsite
##                  diff        lwr         upr     p adj
## shrub-open -0.5018942 -0.9167114 -0.08707699 0.0180668
TukeyHSD(m2, "Nutrient")
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Microsite, arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Nutrient, arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Microsite, Nutrient, arid.gradient
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Plantago.biomass) ~ Microsite * Nutrient * arid.gradient, data = subset(census.arid, Plantago.biomass > 0))
## 
## $Nutrient
##                     diff      lwr        upr     p adj
## Ambient-Added -0.8205043 -1.23339 -0.4076185 0.0001319
## run full model for salvia biomass
m3 <- aov(log(Salvia.biomass)~ Microsite * Nutrient * arid.gradient,  data=subset(census.arid, Salvia.biomass>0))
summary(m3) ## site and year + sitexmicrosite
##                                   Df Sum Sq Mean Sq F value   Pr(>F)    
## Microsite                          1    0.1    0.11   0.054 0.815771    
## Nutrient                           1   24.8   24.84  12.494 0.000465 ***
## arid.gradient                      1   33.7   33.73  16.965 4.79e-05 ***
## Microsite:Nutrient                 1    0.2    0.23   0.115 0.734281    
## Microsite:arid.gradient            1    0.6    0.59   0.299 0.585029    
## Nutrient:arid.gradient             1   20.3   20.28  10.199 0.001538 ** 
## Microsite:Nutrient:arid.gradient   1    0.0    0.00   0.000 0.988929    
## Residuals                        338  672.1    1.99                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m3, "Microsite")
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Microsite, arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Nutrient, arid.gradient
## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## Microsite, Nutrient, arid.gradient
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Salvia.biomass) ~ Microsite * Nutrient * arid.gradient, data = subset(census.arid, Salvia.biomass > 0))
## 
## $Microsite
##                  diff      lwr       upr     p adj
## shrub-open 0.03535378 -0.26289 0.3335975 0.8157715
m3.n <- lm(log(Salvia.biomass)~ arid.gradient,  data=subset(census.arid, Salvia.biomass>0 & Nutrient =="Added"))
m3.a <- lm(log(Salvia.biomass)~ arid.gradient,  data=subset(census.arid, Salvia.biomass>0 & Nutrient =="Ambient"))

## abundance
## run full model for phacelia abundance
m1 <- glm.nb(Phacelia~ Microsite * Nutrient * arid.gradient,  data=census.arid)
summary(m1) ## site and microsite
## 
## Call:
## glm.nb(formula = Phacelia ~ Microsite * Nutrient * arid.gradient, 
##     data = census.arid, init.theta = 0.2308296117, link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.03830  -0.81531  -0.76103   0.01823   3.03660  
## 
## Coefficients:
##                                                Estimate Std. Error z value
## (Intercept)                                  -0.0007756  0.5883363  -0.001
## Micrositeshrub                                1.2705913  0.8115605   1.566
## NutrientAmbient                               0.1077404  0.8333728   0.129
## arid.gradient                                 0.1858087  0.1944620   0.956
## Micrositeshrub:NutrientAmbient                0.1849576  1.1409892   0.162
## Micrositeshrub:arid.gradient                  0.2621592  0.2696848   0.972
## NutrientAmbient:arid.gradient                 0.0421013  0.2759568   0.153
## Micrositeshrub:NutrientAmbient:arid.gradient -0.0505405  0.3789622  -0.133
##                                              Pr(>|z|)
## (Intercept)                                     0.999
## Micrositeshrub                                  0.117
## NutrientAmbient                                 0.897
## arid.gradient                                   0.339
## Micrositeshrub:NutrientAmbient                  0.871
## Micrositeshrub:arid.gradient                    0.331
## NutrientAmbient:arid.gradient                   0.879
## Micrositeshrub:NutrientAmbient:arid.gradient    0.894
## 
## (Dispersion parameter for Negative Binomial(0.2308) family taken to be 1)
## 
##     Null deviance: 585.87  on 839  degrees of freedom
## Residual deviance: 550.93  on 832  degrees of freedom
## AIC: 1941.5
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.2308 
##           Std. Err.:  0.0231 
## 
##  2 x log-likelihood:  -1923.4540
anova(m1, test="Chisq")
## Warning in anova.negbin(m1, test = "Chisq"): tests made without re-
## estimating 'theta'
## Analysis of Deviance Table
## 
## Model: Negative Binomial(0.2308), link: log
## 
## Response: Phacelia
## 
## Terms added sequentially (first to last)
## 
## 
##                                  Df Deviance Resid. Df Resid. Dev
## NULL                                               839     585.87
## Microsite                         1  19.9145       838     565.96
## Nutrient                          1   1.0729       837     564.89
## arid.gradient                     1  11.4485       836     553.44
## Microsite:Nutrient                1   1.0254       835     552.41
## Microsite:arid.gradient           1   1.4563       834     550.96
## Nutrient:arid.gradient            1   0.0064       833     550.95
## Microsite:Nutrient:arid.gradient  1   0.0167       832     550.93
##                                   Pr(>Chi)    
## NULL                                          
## Microsite                        8.098e-06 ***
## Nutrient                         0.3002852    
## arid.gradient                    0.0007155 ***
## Microsite:Nutrient               0.3112354    
## Microsite:arid.gradient          0.2275227    
## Nutrient:arid.gradient           0.9362121    
## Microsite:Nutrient:arid.gradient 0.8971935    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lsmeans(m1, pairwise~"Microsite")
## NOTE: Results may be misleading due to involvement in interactions
## $lsmeans
##  Microsite     lsmean        SE df  asymp.LCL  asymp.UCL
##  open      -0.5548206 0.1205597 NA -0.7911133 -0.3185280
##  shrub      0.1120837 0.1124787 NA -0.1083705  0.3325379
## 
## Results are averaged over the levels of: Nutrient 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast       estimate        SE df z.ratio p.value
##  open - shrub -0.6669043 0.1648821 NA  -4.045  0.0001
## 
## Results are averaged over the levels of: Nutrient 
## Results are given on the log (not the response) scale.
## run full model for plantago biomass
m2 <- glm.nb(Plantago~ Microsite * Nutrient * arid.gradient,  data=census.arid)
summary(m2) ## site and microsite
## 
## Call:
## glm.nb(formula = Plantago ~ Microsite * Nutrient * arid.gradient, 
##     data = census.arid, init.theta = 0.1287536289, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8053  -0.7211  -0.6578  -0.5689   2.4129  
## 
## Coefficients:
##                                                Estimate Std. Error z value
## (Intercept)                                   1.0031515  0.7094622   1.414
## Micrositeshrub                               -2.8208986  1.0387678  -2.716
## NutrientAmbient                              -0.0810403  1.0160835  -0.080
## arid.gradient                                 0.3420451  0.2337751   1.463
## Micrositeshrub:NutrientAmbient                0.5580731  1.4672180   0.380
## Micrositeshrub:arid.gradient                 -0.7058890  0.3368271  -2.096
## NutrientAmbient:arid.gradient                 0.0826214  0.3361349   0.246
## Micrositeshrub:NutrientAmbient:arid.gradient -0.0006522  0.4776644  -0.001
##                                              Pr(>|z|)   
## (Intercept)                                   0.15737   
## Micrositeshrub                                0.00662 **
## NutrientAmbient                               0.93643   
## arid.gradient                                 0.14343   
## Micrositeshrub:NutrientAmbient                0.70368   
## Micrositeshrub:arid.gradient                  0.03611 * 
## NutrientAmbient:arid.gradient                 0.80584   
## Micrositeshrub:NutrientAmbient:arid.gradient  0.99891   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.1288) family taken to be 1)
## 
##     Null deviance: 427.38  on 839  degrees of freedom
## Residual deviance: 414.78  on 832  degrees of freedom
## AIC: 1596.1
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.1288 
##           Std. Err.:  0.0138 
## 
##  2 x log-likelihood:  -1578.1100
anova(m2, test="Chisq")
## Warning in anova.negbin(m2, test = "Chisq"): tests made without re-
## estimating 'theta'
## Analysis of Deviance Table
## 
## Model: Negative Binomial(0.1288), link: log
## 
## Response: Plantago
## 
## Terms added sequentially (first to last)
## 
## 
##                                  Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                                               839     427.38         
## Microsite                         1   5.5584       838     421.82  0.01839
## Nutrient                          1   0.0296       837     421.79  0.86334
## arid.gradient                     1   0.5610       836     421.23  0.45386
## Microsite:Nutrient                1   1.6621       835     419.56  0.19732
## Microsite:arid.gradient           1   4.7218       834     414.84  0.02978
## Nutrient:arid.gradient            1   0.0649       833     414.78  0.79889
## Microsite:Nutrient:arid.gradient  1   0.0000       832     414.78  0.99859
##                                   
## NULL                              
## Microsite                        *
## Nutrient                          
## arid.gradient                     
## Microsite:Nutrient                
## Microsite:arid.gradient          *
## Nutrient:arid.gradient            
## Microsite:Nutrient:arid.gradient  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lsmeans(m2, pairwise~"Microsite")
## NOTE: Results may be misleading due to involvement in interactions
## $lsmeans
##  Microsite     lsmean        SE df  asymp.LCL  asymp.UCL
##  open      -0.1639688 0.1466943 NA -0.4514843  0.1235467
##  shrub     -0.6304166 0.1522519 NA -0.9288248 -0.3320084
## 
## Results are averaged over the levels of: Nutrient 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast      estimate        SE df z.ratio p.value
##  open - shrub 0.4664478 0.2114234 NA   2.206  0.0274
## 
## Results are averaged over the levels of: Nutrient 
## Results are given on the log (not the response) scale.
## run full model for salvia biomass
m3 <- glm.nb(Salvia~ Microsite * Nutrient * arid.gradient,  data=census.arid)
summary(m3) ## site and microsite
## 
## Call:
## glm.nb(formula = Salvia ~ Microsite * Nutrient * arid.gradient, 
##     data = census.arid, init.theta = 0.3061627656, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3585  -1.1226  -0.8020   0.2021   2.2024  
## 
## Coefficients:
##                                              Estimate Std. Error z value
## (Intercept)                                   2.39712    0.46918   5.109
## Micrositeshrub                               -0.02997    0.67161  -0.045
## NutrientAmbient                               0.61208    0.67236   0.910
## arid.gradient                                 0.57763    0.15655   3.690
## Micrositeshrub:NutrientAmbient                0.16062    0.94881   0.169
## Micrositeshrub:arid.gradient                  0.07096    0.22516   0.315
## NutrientAmbient:arid.gradient                 0.24763    0.22629   1.094
## Micrositeshrub:NutrientAmbient:arid.gradient -0.14470    0.31899  -0.454
##                                              Pr(>|z|)    
## (Intercept)                                  3.24e-07 ***
## Micrositeshrub                               0.964406    
## NutrientAmbient                              0.362646    
## arid.gradient                                0.000225 ***
## Micrositeshrub:NutrientAmbient               0.865571    
## Micrositeshrub:arid.gradient                 0.752644    
## NutrientAmbient:arid.gradient                0.273802    
## Micrositeshrub:NutrientAmbient:arid.gradient 0.650108    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.3062) family taken to be 1)
## 
##     Null deviance: 785.98  on 839  degrees of freedom
## Residual deviance: 721.96  on 832  degrees of freedom
## AIC: 2994
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.3062 
##           Std. Err.:  0.0230 
## 
##  2 x log-likelihood:  -2976.0340
anova(m3, test="Chisq")
## Warning in anova.negbin(m3, test = "Chisq"): tests made without re-
## estimating 'theta'
## Analysis of Deviance Table
## 
## Model: Negative Binomial(0.3062), link: log
## 
## Response: Salvia
## 
## Terms added sequentially (first to last)
## 
## 
##                                  Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                                               839     785.98         
## Microsite                         1    0.258       838     785.72  0.61139
## Nutrient                          1    3.609       837     782.11  0.05747
## arid.gradient                     1   54.596       836     727.51 1.48e-13
## Microsite:Nutrient                1    4.529       835     722.98  0.03332
## Microsite:arid.gradient           1    0.005       834     722.98  0.94599
## Nutrient:arid.gradient            1    0.870       833     722.11  0.35089
## Microsite:Nutrient:arid.gradient  1    0.145       832     721.96  0.70327
##                                     
## NULL                                
## Microsite                           
## Nutrient                         .  
## arid.gradient                    ***
## Microsite:Nutrient               *  
## Microsite:arid.gradient             
## Nutrient:arid.gradient              
## Microsite:Nutrient:arid.gradient    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lsmeans(m3, pairwise~Microsite*Nutrient)
## NOTE: Results may be misleading due to involvement in interactions
## $lsmeans
##  Microsite Nutrient    lsmean        SE df asymp.LCL asymp.UCL
##  open      Added    0.6995864 0.1351555 NA 0.4346865 0.9644863
##  shrub     Added    0.4610742 0.1382722 NA 0.1900656 0.7320828
##  open      Ambient  0.5839160 0.1380287 NA 0.3133848 0.8544472
##  shrub     Ambient  0.9312613 0.1338667 NA 0.6688874 1.1936351
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                       estimate        SE df z.ratio p.value
##  open,Added - shrub,Added      0.2385122 0.1933552 NA   1.234  0.6054
##  open,Added - open,Ambient     0.1156704 0.1931811 NA   0.599  0.9325
##  open,Added - shrub,Ambient   -0.2316749 0.1902296 NA  -1.218  0.6154
##  shrub,Added - open,Ambient   -0.1228418 0.1953743 NA  -0.629  0.9228
##  shrub,Added - shrub,Ambient  -0.4701871 0.1924565 NA  -2.443  0.0692
##  open,Ambient - shrub,Ambient -0.3473453 0.1922816 NA  -1.806  0.2702
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates
## treatment plots
census.long <- gather(census.arid, species, biomass, 12:14)
census.long[,"species"] <- ifelse(census.long[,"species"] =="Phacelia.biomass", "P. tanacetifolia",
                                  ifelse(census.long[,"species"]=="Plantago.biomass","P. insularis","S. columbariae"))
census.long[,"species"] <- as.factor(census.long$species)

## microsite average
census.micro <- census.long %>% filter(biomass>0) %>%  group_by(Microsite, species) %>%  summarize(Biomass=mean(biomass), error=se(biomass))

ggplot(census.micro, aes(x=species, y=Biomass, fill=Microsite)) +   geom_bar(position=position_dodge(), stat="identity", colour="black") + scale_fill_brewer()+  geom_errorbar(aes(ymin=Biomass-error, ymax=Biomass+error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+ theme_Publication() + theme(axis.text.x = element_text(size=12, face="italic"))
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## nutrient average
census.nut <- census.long %>% filter(biomass>0) %>%  group_by(Nutrient, species, Year) %>%  summarize(Biomass=mean(biomass), error=se(biomass))

ggplot(census.nut, aes(x=species, y=Biomass, fill=Nutrient)) +   geom_bar(position=position_dodge(), stat="identity",colour="black") + scale_fill_brewer()+  geom_errorbar(aes(ymin=Biomass-error, ymax=Biomass+error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+ theme_Publication()+ theme(axis.text.x = element_text(size=12, face="italic")) +facet_grid(~Year)
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## aridity average

ggplot(census.long, aes(x=arid.gradient, y=biomass, fill=species, colour=species)) + theme_Publication()+   stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=biomass), color="#009E73", fill="#009E7340", data=subset(census.long, species=="Phacelia.biomass" & biomass>0)) + stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=biomass), color="#E69F00",fill="#E69F0040", data=subset(census.long, species=="Plantago.biomass" & biomass>0))+ stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=biomass), color="#56B4E9", fill="#56B4E940", data=subset(census.long, species=="Salvia.biomass" & biomass>0))
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

### need to separate zeros using gamma-hurdle model
##http://seananderson.ca/2014/05/18/gamma-hurdle.html


## logistic first
census.long[,"presence"] <- ifelse(census.long$biomass>0,1,0)

m1.l <- glmer(presence~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year),data=subset(census.long, species=="P. tanacetifolia"), family=binomial)
car::Anova(m1.l)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: presence
##                                             Chisq Df Pr(>Chisq)    
## Microsite                                 28.8440  1  7.845e-08 ***
## poly(arid.gradient, 2)                    13.7669  2   0.001025 ** 
## Nutrient                                   0.0065  1   0.935820    
## Microsite:poly(arid.gradient, 2)           1.6297  2   0.442697    
## Microsite:Nutrient                         1.9145  1   0.166465    
## poly(arid.gradient, 2):Nutrient            0.5605  2   0.755599    
## Microsite:poly(arid.gradient, 2):Nutrient  1.5360  2   0.463938    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- lmer(log(biomass)~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year),data=subset(census.long, species=="P. tanacetifolia" & presence==1))
shapiro.test(resid(m1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m1)
## W = 0.99429, p-value = 0.5244
anova(m1)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                                           Sum Sq Mean Sq NumDF  DenDF
## Microsite                                 10.828  10.828     1 220.60
## poly(arid.gradient, 2)                    96.805  48.403     2 193.15
## Nutrient                                   7.785   7.785     1 220.10
## Microsite:poly(arid.gradient, 2)           3.723   1.862     2 220.61
## Microsite:Nutrient                         0.068   0.068     1 220.08
## poly(arid.gradient, 2):Nutrient            4.093   2.046     2 220.02
## Microsite:poly(arid.gradient, 2):Nutrient  2.078   1.039     2 220.06
##                                           F.value   Pr(>F)    
## Microsite                                  5.6947  0.01786 *  
## poly(arid.gradient, 2)                    25.4555 1.54e-10 ***
## Nutrient                                   4.0941  0.04424 *  
## Microsite:poly(arid.gradient, 2)           0.9790  0.37731    
## Microsite:Nutrient                         0.0357  0.85037    
## poly(arid.gradient, 2):Nutrient            1.0762  0.34269    
## Microsite:poly(arid.gradient, 2):Nutrient  0.5464  0.57983    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m1) ## R squared values 
##       R2m       R2c 
## 0.2463743 0.3096945
m2.l <- glmer(presence~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year),data=subset(census.long, species=="P. insularis"), family=binomial)
car::Anova(m2.l)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: presence
##                                             Chisq Df Pr(>Chisq)    
## Microsite                                  1.1386  1   0.285954    
## poly(arid.gradient, 2)                    24.7326  2   4.26e-06 ***
## Nutrient                                   0.3030  1   0.582001    
## Microsite:poly(arid.gradient, 2)           9.8407  2   0.007296 ** 
## Microsite:Nutrient                         3.7231  1   0.053664 .  
## poly(arid.gradient, 2):Nutrient            9.3310  2   0.009415 ** 
## Microsite:poly(arid.gradient, 2):Nutrient  1.8086  2   0.404829    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2 <- lmer(log(biomass)~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year), data=subset(census.long, species=="P. insularis" & presence==1))
shapiro.test(resid(m2))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2)
## W = 0.9915, p-value = 0.49
anova(m2)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                                            Sum Sq Mean Sq NumDF  DenDF
## Microsite                                  0.4051  0.4051     1 141.10
## poly(arid.gradient, 2)                    13.1478  6.5739     2 141.72
## Nutrient                                   7.8434  7.8434     1 141.28
## Microsite:poly(arid.gradient, 2)           0.4061  0.2031     2 141.30
## Microsite:Nutrient                         0.4193  0.4193     1 141.23
## poly(arid.gradient, 2):Nutrient            1.6505  0.8252     2 141.24
## Microsite:poly(arid.gradient, 2):Nutrient  0.3606  0.1803     2 141.17
##                                           F.value   Pr(>F)   
## Microsite                                  0.3818 0.537661   
## poly(arid.gradient, 2)                     6.1948 0.002635 **
## Nutrient                                   7.3911 0.007377 **
## Microsite:poly(arid.gradient, 2)           0.1914 0.826049   
## Microsite:Nutrient                         0.3951 0.530642   
## poly(arid.gradient, 2):Nutrient            0.7777 0.461438   
## Microsite:poly(arid.gradient, 2):Nutrient  0.1699 0.843906   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values 
##       R2m       R2c 
## 0.1291956 0.6996902
m3.l <- glmer(presence~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year),data=subset(census.long, species=="S. columbariae"), family=binomial)
car::Anova(m3.l)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: presence
##                                             Chisq Df Pr(>Chisq)    
## Microsite                                  0.0208  1    0.88536    
## poly(arid.gradient, 2)                    72.4822  2    < 2e-16 ***
## Nutrient                                   0.9592  1    0.32739    
## Microsite:poly(arid.gradient, 2)           2.2173  2    0.33000    
## Microsite:Nutrient                         3.8296  1    0.05036 .  
## poly(arid.gradient, 2):Nutrient            2.1024  2    0.34951    
## Microsite:poly(arid.gradient, 2):Nutrient  0.3469  2    0.84075    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 <- lmer(log(Salvia.biomass)~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year), data=subset(census.arid, Salvia.biomass>0))
shapiro.test(resid(m3))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m3)
## W = 0.99339, p-value = 0.1337
anova(m3)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                                           Sum Sq Mean Sq NumDF  DenDF
## Microsite                                  1.519  1.5187     1 333.08
## poly(arid.gradient, 2)                    52.112 26.0561     2 333.65
## Nutrient                                  15.470 15.4701     1 333.09
## Microsite:poly(arid.gradient, 2)           1.015  0.5076     2 333.15
## Microsite:Nutrient                         0.244  0.2444     1 333.03
## poly(arid.gradient, 2):Nutrient           17.550  8.7750     2 333.04
## Microsite:poly(arid.gradient, 2):Nutrient  4.472  2.2362     2 333.06
##                                           F.value    Pr(>F)    
## Microsite                                  0.9019  0.342966    
## poly(arid.gradient, 2)                    15.4738 3.745e-07 ***
## Nutrient                                   9.1871  0.002628 ** 
## Microsite:poly(arid.gradient, 2)           0.3014  0.739951    
## Microsite:Nutrient                         0.1451  0.703473    
## poly(arid.gradient, 2):Nutrient            5.2112  0.005909 ** 
## Microsite:poly(arid.gradient, 2):Nutrient  1.3280  0.266399    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m3) ## R squared values 
##       R2m       R2c 
## 0.1322110 0.2743503
## set up prediction dataset
library(effects)
## Loading required package: carData
## Use the command
##     lattice::trellis.par.set(effectsTheme())
##   to customize lattice options for effects plots.
## See ?efffectTheme for details.
epha <- as.data.frame(Effect(c("arid.gradient"),m1.l))
epla <- as.data.frame(Effect(c("arid.gradient"),m2.l))
esal <- as.data.frame(Effect(c("arid.gradient"),m3.l))
ee <- rbind(epha,epla, esal)
ee[,"species"] <- rep(unique(census.long$species),each=5)


ggplot(data=ee, aes(arid.gradient,fit,color=species))+geom_line(lwd=1.4) +xlab("aridity gradient") + ylab("probability of occurrence") + theme_Publication() +  geom_ribbon(aes(ymin=lower, ymax=upper, fill=species),  alpha=0.5)
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## Calculate effects size estimates for regressions
library(effects)
epha <- as.data.frame(Effect(c("arid.gradient"),m1))
epla <- as.data.frame(Effect(c("arid.gradient"),m2))
esal <- as.data.frame(Effect(c("arid.gradient"),m3))
ee <- rbind(epha,epla, esal)
ee[,"species"] <- rep(unique(census.long$species),each=5)
#ee[,2:5] <- exp(ee[,2:5])

## add colours for species
color.add<- subset(census.long, biomass>0)
color.add$color <- ifelse(color.add$species=="P. tanacetifolia", "#009E73",ifelse(color.add$species=="P. insularis", "#D55E00","#0072B2"))

ggplot()+ geom_jitter(aes(x=arid.gradient, y=log(biomass), group=species), color=color.add$color,alpha=0.4, width = 0.05, data=subset(census.long, biomass>0)) +geom_line(data=ee, aes(arid.gradient,fit,color=species),lwd=1.4) +xlab("aridity gradient") + ylab("log-transformed biomass") + theme_Publication()
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

# 
# ## to check for overdispersion
# overdisp_fun <- function(model) {
#   ## number of variance parameters in 
#   ##   an n-by-n variance-covariance matrix
#   vpars <- function(m) {
#     nrow(m)*(nrow(m)+1)/2
#   }
#   model.df <- sum(sapply(VarCorr(model),vpars))+length(fixef(model))
#   (rdf <- nrow(model@frame)-model.df)
#   rp <- residuals(model)
#   Pearson.chisq <- sum(rp^2)
#   prat <- Pearson.chisq/rdf
#   pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE,log.p=TRUE)
#   c(chisq=Pearson.chisq,ratio=prat,p=exp(pval))
# }


## related packages
# library(car)
# library(lme4)


## https://www.rdocumentation.org/packages/lme4/versions/1.1-15/topics/pvalues 
## https://arxiv.org/pdf/1406.5823.pdf
## Calculating Pvalues from glmer

# ##Provides parametric bootstrapping methods for model comparisons
# library(afex)
# mixed(Phacelia~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year), data=census.arid, family=poisson, control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)), method="LRT")


# ## accounting for overdispersion in GLMM
# ## Harrison, X. A. (2014). Using observation-level random effects to model overdispersion in count data in ecology and evolution. PeerJ, 2, e616.
# m4 <- glmer(Phacelia~Microsite * arid.gradient * Nutrient + (1|Year), data=census.arid, family=poisson)
# #car::Anova(m4, test.statistics="LR", type=c("III")) 
# overdisp_fun(m4)  ## over dispersion?
# 
# 
# ## re-run with observational-level random effects
# disp<-seq(nrow(census.arid))
# m4 <- glmer(Phacelia~Microsite * arid.gradient * Nutrient + (1|Year) + (1|disp), data=census.arid, family=poisson)
# car::Anova(m4, test.statistics="LR", type=c("III")) 
# 
# m4 <- glmer.nb(Phacelia~Microsite * arid.gradient * Nutrient + (1|Year), data=census.arid)
# car::Anova(m4, test.statistics="LR", type=c("III")) 
# 
# m5 <- glmer(Plantago~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year), data=census.arid, family=poisson)
# overdisp_fun(m5)  ## over dispersion?
# 
# ## re-run with observational-level random effects
# m5 <- glmer(Plantago~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year)+ (1|disp), data=census.arid, family=poisson)
# car::Anova(m5, test.statistics="LR", type=c("III"))
# 
# m6 <- glmer(Salvia~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year), data=census.arid, family=poisson)
# overdisp_fun(m6)  ## over dispersion?
# 
# 
# ## re-run with observational-level random effects
# m6 <- glmer(Salvia~Microsite * poly(arid.gradient,2) * Nutrient + (1|Year)+ (1|disp), data=census.arid, family=poisson)
# car::Anova(m6, test.statistics="LR", type=c("III"))

Community analysis

## Regressions To Separate Phylogenetic Attraction And Repulsion

## sppregs(community[,13:53], )

### Network analysis
library(igraph)
library(statnet)


## reload native vs non-native species
status <- read.csv("Data//ERG.specieslist.csv")

shrub.comm <- subset(community, Microsite=="shrub", select=13:53)
shrub.comm <- shrub.comm[,colSums(shrub.comm)!=0]
open.comm <- subset(community, Microsite=="open", select=13:53)
open.comm <- open.comm[,colSums(open.comm)!=0]

## Shrub
## correlation species network
network <- data.frame(cor(shrub.comm))

network[,"To"] <- colnames(network)

test <- gather(network, From, correlation, 1:(ncol(network)-1))
test <- subset(test, correlation>0.10 & correlation != 1 | correlation < -0.10 )

n.shrub <- graph_from_data_frame(test[,1:2], directed=TRUE)

## plot network circle
lo <- layout_in_circle(n.shrub)
plot(n.shrub, layout=lo)
transitivity(n.shrub, type = "local") 
transitivity(n.shrub, type = "global") 


## colour graphs
shapes <- V(n.shrub)$name %in% (subset(status, status=="non.native", select=1))[,1] +1
plot(n.shrub, layout=layout_on_grid(n.shrub),vertex.shape=c("circle", "square")[shapes],  vertex.color=c("tomato2", "royalblue")[shapes])

## open
## correlation species network
network <- data.frame(cor(open.comm))

network[,"To"] <- colnames(network)

test <- gather(network, From, correlation, 1:(ncol(network)-1))
test <- subset(test, correlation>0.10 & correlation != 1 | correlation < -0.10 )

n.open <- graph_from_data_frame(test[,1:2], directed=TRUE)

## plot network circle
lo <- layout_in_circle(n.open)
plot(n.open, layout=lo)

## colour graphs
shapes <- V(n.open)$name %in% (subset(status, status=="non.native", select=1))[,1] +1
plot(n.open, layout=layout_on_grid(n.open),vertex.shape=c("circle", "square")[shapes],  vertex.color=c("tomato2", "royalblue")[shapes])

## statistics
transitivity(n.open, type = "local") 
transitivity(n.open, type = "global") 

indicator species analysis

library(indicspecies)

## indicator species analysis
indval <- multipatt(community[,13:53], community$Microsite, control=how(nperm=999))
summary(indval)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: IndVal.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 41
##  Selected number of species: 14 
##  Number of species associated to 1 group: 14 
## 
##  List of species associated to each combination: 
## 
##  Group open  #sps.  8 
##                   stat p.value    
## E.cicutarium     0.558   0.001 ***
## L.gracilis       0.427   0.001 ***
## A.wrangelianus   0.201   0.001 ***
## C.brevipes       0.171   0.029 *  
## A.lentiginosus   0.149   0.030 *  
## C.rigida         0.149   0.016 *  
## A.grandiflora    0.129   0.014 *  
## Boraginaceae.sp. 0.120   0.018 *  
## 
##  Group shrub  #sps.  6 
##                  stat p.value    
## C.fremontii     0.379   0.001 ***
## M.glabrata      0.277   0.001 ***
## P.tanacetifolia 0.267   0.001 ***
## A.tessellata    0.199   0.002 ** 
## B.nigra         0.179   0.012 *  
## E.nitidum       0.172   0.008 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
site.names <- as.character(unique(community$Site))
community.site <- subset(community, Site==site.names[7])



## species accumilation curves
accum1 <- specaccum(subset(community.site, Microsite=="shrub", select=13:53), method="random",  permutations=1000)
accum2 <- specaccum(subset(community.site, Microsite=="open", select=13:53), method="random",  permutations=1000)


par(mar=c(4.5,4.5,0.5,0.5))
plot(accum2,  col="#E69F00",ci.col="#E69F0090", xlab="Sampling effort", cex.axis=1.5, cex.lab=1.8, ylab="Species richness", lwd=2, ylim=c(0,16))
plot(accum1,add=T , col="#56B4E9",   ci.col="#56B4E990", lwd=2)

micros <- community$Microsite[1:120]

indval <- multipatt(community.site[,13:53], micros, control=how(nperm=999))
summary(indval)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: IndVal.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 41
##  Selected number of species: 4 
##  Number of species associated to 1 group: 4 
## 
##  List of species associated to each combination: 
## 
##  Group open  #sps.  1 
##                   stat p.value  
## Boraginaceae.sp. 0.316   0.025 *
## 
##  Group shrub  #sps.  3 
##               stat p.value    
## A.tessellata 0.531   0.001 ***
## B.nigra      0.475   0.012 *  
## M.affinis    0.343   0.046 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cluster analysis

library("cluster")
library("factoextra")
library("magrittr")

comm.sum <- community %>%  group_by(Microsite, Site) %>% summarize_if(is.numeric, funs(sum))
comm.sum <- data.frame(comm.sum)
rownames(comm.sum) <- paste(comm.sum$Microsite,comm.sum$Site)

## coreelation matrix

### Distance measures
res.dist <- get_dist(comm.sum[,13:53], stand = TRUE, method = "euclidean")
fviz_dist(res.dist,    gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))

## dendrogram

## transform data
comm.trans <- decostand(comm.sum[,13:53], "hell")

## clean data for ordination
## see distribution of spp
boxplot(comm.trans, xaxt="n")
labs <- colnames(comm.trans)
text(cex=0.8, x=1:41-1, y=-0.12, labs, xpd=TRUE, srt=45)

## remove spp with only one instancve
comm.trans <- comm.trans[,!colSums(comm.trans)==apply(comm.trans, 2, max)]


## replace outliers with mean
avg.max <- function(x) {
  y =  max(x)
  avg = mean(x)
  x[x==y] <- avg
  return(x)
}
  
## dot chart to find outliers (remove or Avg)
#dotchart(comm.trans$E.cicutarium) ## Erodium
comm.trans[,"E.cicutarium"] <- avg.max(comm.trans$E.cicutarium) ## removed outlier
#dotchart(comm.trans$A.wrangelianus) ## Chilean trefoil
comm.trans[,"A.wrangelianus"] <- NULL ## removed entire species
#dotchart(comm.trans$A.lentiginosus) ##
comm.trans[,"A.lentiginosus"] <- avg.max(comm.trans$A.lentiginosus) ## removed outlier
#dotchart(comm.trans$H.vulgare) ## Hordeum
comm.trans[,"H.vulgare"] <- NULL ## removed entire species
#dotchart(comm.trans$Pectocarya.spp) ## 
comm.trans[,"Pectocarya.spp"] <- avg.max(comm.trans$Pectocarya.spp) ## removed outlier
#dotchart(comm.trans$L.arizonicus) ## Lupin 
comm.trans[,"L.arizonicus"] <- NULL ## removed entire species
#dotchart(comm.trans$M.bellioides) ##
comm.trans[,"M.bellioides"] <- NULL ## removed entire species
#dotchart(comm.trans$B.nigra) ## Mustard
comm.trans[,"B.nigra"] <- NULL ## removed entire species
#dotchart(comm.trans$N.demissum) ## Sand Matt
comm.trans[,"N.demissum"] <- avg.max(comm.trans$N.demissum) ## removed outlier
#dotchart(comm.trans$L.gracilis) ## 
comm.trans[,"L.gracilis"] <- avg.max(comm.trans$L.gracilis) ## removed outlier
#dotchart(comm.trans$A.grandiflora) ##  Amsinckia
comm.trans[,"A.grandiflora"] <- NULL ## removed entire species
#dotchart(comm.trans$B.diandrus) ##  Bromus diandrus
comm.trans[,"B.diandrus"] <- NULL ## removed entire species
#dotchart(comm.trans$M.affinis) ##  
comm.trans[,"M.affinis"] <- NULL ## removed entire species
#dotchart(comm.trans$P.crenulata) ##  
comm.trans[,"P.crenulata"] <- NULL ## removed entire species
#dotchart(comm.trans$Erinoginum.spp) # buckwheat
comm.trans[,"Erinoginum.spp"] <- avg.max(comm.trans$Erinoginum.spp) ## removed outlier
#dotchart(comm.trans$C.lasiophyllus) # 
comm.trans[,"C.lasiophyllus"] <- NULL ## removed entire species



library(usdm)
library(corrplot)

## Check for collinearity
vifcor(comm.trans)
## 4 variables from the 25 input variables have collinearity problem: 
##  
## E.wallacei P.insularis E.glyptosperma L.decorus 
## 
## After excluding the collinear variables, the linear correlation coefficients ranges between: 
## min correlation ( L.gracilis ~ S.barbatus ):  -0.003266766 
## max correlation ( C.fremontii ~ M.glabrata ):  0.8450794 
## 
## ---------- VIFs of the remained variables -------- 
##          Variables VIF
## 1         B.rubens Inf
## 2     E.cicutarium Inf
## 3       S.barbatus Inf
## 4      C.barbigera Inf
## 5    C.claviformis Inf
## 6   A.lentiginosus Inf
## 7         C.rigida Inf
## 8  P.tanacetifolia Inf
## 9   Pectocarya.spp Inf
## 10     D.capitatum Inf
## 11      C.brevipes Inf
## 12      M.glabrata Inf
## 13     C.fremontii Inf
## 14    C.intermedia Inf
## 15      N.demissum Inf
## 16      L.gracilis Inf
## 17    A.tessellata Inf
## 18       P.secunda Inf
## 19  L.sparsiflorus Inf
## 20  Erinoginum.spp Inf
## 21       E.nitidum Inf
## drop Claviformis because of collinear problems in P.insularis
comm.trans[,"C.claviformis"] <- NULL ## removed entire species
comm.trans[,"L.decorus"] <- NULL ## removed entire species
comm.trans[,"E.glyptosperma"] <- NULL ## removed entire species
comm.trans[,"E.wallacei"] <- NULL ## removed entire species

cor.dat <- cor(comm.trans[,1:20])
corrplot(cor.dat, method="number")

## calculate distances
dis <- vegdist(comm.trans, method="bray")


## cluster analysis
clus <- hclust(dis, "ward.D2")

## cluster colours
colz <- c("#3e4444", "#86af49", "#405d27","#c1946a")

fviz_dend(clus, k = 4, # Cut in four groups
          cex = 0.7, cex.axis=2, # label size
          k_colors = colz,
          color_labels_by_k = TRUE, # color labels by groups
          rect = TRUE # Add rectangle around groups
          ) 

grp <- cutree(clus, 4)

## Interpretation of classes

## get soil moisture from site/microsite
smc.avg <- smc.end %>% group_by(Site, Microsite) %>% summarize(mean.smc=mean(swc.avg, na.rm=T))
smc.avg <- data.frame(smc.avg)
smc.avg[,"grp"] <- as.factor(c(3,3,2,1,4,4,1,1,4,4,3,3,2,2))

boxplot(mean.smc ~ grp, data=smc.avg)

## CA or PCA
dca1 <- decorana(comm.trans) ## length of gradient >2 & determine relative differences in community composition 


## conduct ordination
ord <- cca(comm.trans)
summary(ord)
## 
## Call:
## cca(X = comm.trans) 
## 
## Partitioning of mean squared contingency coefficient:
##               Inertia Proportion
## Total            1.53          1
## Unconstrained    1.53          1
## 
## Eigenvalues, and their contribution to the mean squared contingency coefficient 
## 
## Importance of components:
##                          CA1    CA2    CA3     CA4     CA5     CA6     CA7
## Eigenvalue            0.6615 0.2873 0.1656 0.10566 0.08765 0.07595 0.05451
## Proportion Explained  0.4324 0.1878 0.1082 0.06907 0.05729 0.04964 0.03563
## Cumulative Proportion 0.4324 0.6202 0.7284 0.79749 0.85478 0.90443 0.94006
##                           CA8     CA9    CA10     CA11    CA12     CA13
## Eigenvalue            0.03070 0.02546 0.01821 0.008974 0.00457 0.003798
## Proportion Explained  0.02006 0.01664 0.01190 0.005870 0.00299 0.002480
## Cumulative Proportion 0.96012 0.97676 0.98866 0.994530 0.99752 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##                     CA1      CA2      CA3      CA4      CA5       CA6
## B.rubens        -1.5272  0.37620  0.03693 -0.17971  0.26195 -0.007680
## E.cicutarium    -0.5874 -0.19213 -0.20338 -0.20206 -0.36108  0.524976
## S.barbatus       0.1974 -0.64533  0.21879  0.05816  0.09674 -0.053334
## C.barbigera      0.4900 -0.11013  0.11027  0.09640  0.08716  0.152044
## A.lentiginosus   0.2199 -0.28948 -0.59265 -0.08500  0.77156 -0.296433
## C.rigida         0.4907 -0.59656 -0.53938  0.07713  0.37931 -0.042119
## P.tanacetifolia -1.2902  0.91318 -0.45764  2.53829  0.03756 -0.152588
## Pectocarya.spp   0.5495 -0.11384 -0.59153 -0.11011  0.04276  0.166931
## D.capitatum     -0.7173 -0.38053 -0.34301  0.92059 -1.47031  1.938948
## C.brevipes       0.4996 -0.08694 -1.06985 -0.06884 -0.31678  0.438663
## P.insularis      0.8031  0.56576 -0.62483 -0.12389  0.05743 -0.096057
## M.glabrata       0.7884  0.84798  0.94548 -0.09754 -0.18919  0.009159
## C.fremontii      0.8264  0.71519  0.19154 -0.07910 -0.11619  0.016657
## C.intermedia     0.3301 -0.17327 -0.89534 -0.13610 -0.03923 -0.580017
## N.demissum       0.7716  0.44620 -0.72001  0.09033  0.21708 -0.241934
## L.gracilis      -0.9869 -0.31033  0.07421  0.35879 -1.02622 -0.476879
## A.tessellata    -1.0035 -0.36194  0.26967 -0.84774 -1.13864 -1.238336
## P.secunda       -1.5685  0.52156  0.01749  0.80285  0.12080 -0.361127
## L.sparsiflorus   0.7841  1.12673  0.93614  0.24022  0.04660  0.069996
## Erinoginum.spp   0.5728 -1.00733  0.32438  0.36758  1.12435 -0.440374
## E.nitidum        0.8450  0.85547 -0.32831 -0.02364 -0.06814 -0.220852
## 
## 
## Site scores (weighted averages of species scores)
## 
##                           CA1     CA2      CA3      CA4      CA5       CA6
## open Barstow           0.3623 -1.3253  0.59499  0.23046  0.56897  0.613534
## open Cuyama           -0.2820 -0.9387 -0.39810  0.01254 -2.11537  2.774795
## open HeartofMojave     0.7431  0.2335 -2.05624 -0.52149  0.05849  0.033486
## open PanocheHills     -1.6700  0.4724  0.08096 -1.23638  1.43764  0.474226
## open SheepholeValley   0.9186  0.9037  0.77143  0.17618  0.30513  0.116582
## open Tecopa            0.5210 -1.3909 -0.13250  0.41626  1.33826 -0.340390
## open TejonRanch       -0.6455 -0.9594  0.23131 -0.09915 -1.38186 -1.124251
## shrub Barstow          0.5555 -0.5840  0.87451  0.02079  0.01283  0.724468
## shrub Cuyama          -1.7597  0.9562 -0.21109  3.09525  0.07451 -0.062796
## shrub HeartofMojave    0.7511  0.7261 -1.52964  0.11670 -0.12309 -0.542992
## shrub PanocheHills    -2.0371  0.9309 -0.05454 -1.74126  1.62903  1.240252
## shrub SheepholeValley  1.0489  1.6807  1.57641 -0.25781 -0.55721 -0.007623
## shrub Tecopa           0.4692 -1.4952  0.52997  0.40750  1.29037 -0.766954
## shrub TejonRanch      -1.1319 -0.1862  0.35145 -1.12008 -1.13771 -1.730813
## calculate priority
spp.priority <- colSums(comm.trans)

## plot ordination
par(mar=c(4.5,4.5,0.5,0.5))
plot(ord, type="n")
ordiellipse(ord, grp, lty = 2, col = "grey80", draw="polygon", alpha=150)
orditorp(ord, display = "species", cex = 0.7, col = "darkorange3", priority=spp.priority, air=0.5)
orditorp(ord, display = "sites", cex = 0.7, col = "darkslateblue", air=0.1)


##collect environmental variables for site
nutrients <- read.csv("Data/ERG.soilnutrients.csv")
nutrients.mean <- nutrients %>% group_by(microsite, Site) %>% summarise_all(funs(mean))
nutrients.vars <- data.frame(nutrients.mean)


## summarize daily means across sites
means <- HOBO.data %>% group_by(Microsite, Site) %>% summarize(Temp.=mean(Temp, na.rm=T),rh=mean(RH, na.rm=T),temp.se=se(Temp),rh.se=se(RH))
means <- data.frame(means)

## soil moisture
smc.avg <- smc.early %>% group_by(Microsite, Site) %>% summarize(SMC=mean(swc.avg, na.rm=T))

envs <- data.frame(swc=smc.avg[,"SMC"],nutrients.vars[,c("N","P","K")], means[,c("Temp.","rh")])

## Check for collinearity
cor(envs)
##              SMC           N           P           K      Temp.         rh
## SMC    1.0000000 -0.31692385  0.44854712  0.81694831 -0.6306962  0.9176776
## N     -0.3169238  1.00000000  0.05879919 -0.05275422  0.4188434 -0.4374109
## P      0.4485471  0.05879919  1.00000000  0.64885814 -0.5264257  0.5353719
## K      0.8169483 -0.05275422  0.64885814  1.00000000 -0.5705341  0.7740852
## Temp. -0.6306962  0.41884336 -0.52642569 -0.57053407  1.0000000 -0.8625519
## rh     0.9176776 -0.43741089  0.53537193  0.77408521 -0.8625519  1.0000000
envs[,"K"] <- NULL ## drop potassium for being correlated
envs[,"rh"] <- NULL ## drop humidity for being correlated

ord.env <- envfit(ord, envs)
plot(ord.env)

# ordicluster(ord, clus, col="blue")

# ## conduct ordination with gradient analysis
# ord <- cca(comm.trans~SMC+N+P+Temp.,data=envs)
# (R2adj <- RsquareAdj(ord)$adj.r.squared) ## 30.0% of variation explained
# 
# ord <- cca(comm.trans~.,data=envs)
# ord.0 <- cca(comm.trans~1,data=envs)
# m <- ordistep(ord.0, scope=formula(ord))
# 
# 
# ## plot ordination
# par(mar=c(4.5,4.5,0.5,0.5))
# plot(ord, type="n") # display="bp", ylim=c(-2,2),xlim=c(-2,2))
# ordiellipse(ord, grp, lty = 2, col = colz, draw="polygon", alpha=150)
# orditorp(ord, display = "species", cex = 0.7, col = "darkred", priority=spp.priority, air=0.5)
# orditorp(ord, display = "sites", cex = 0.7, col = "black", air=0.1)

Ordination including aridity and microsite (pCCA)

Partial Constrained Correspondence Analysis for community data

Mixed model resuilts

community.arid <- merge(community, gams.data, by=c("Year","Site"))


### Ambient
mean.year <- community.arid %>% group_by(Microsite, gams, Year) %>%  summarize(abd=mean(Abundance),rich=mean(Richness), bio=mean(Biomass))

## combine other dataframes with gams
mpd.gams <- merge(mpd.arid, gams.data, by=c("Year","Site"))
mean.status2 <- merge(mean.spp,gams.data, by=c("Year","Site"))
mean.status2[,"native.ratio"] <- mean.status2$native/(mean.status2$native+mean.status2$non.native)
mean.status2[is.na(mean.status2)] <- 0

m1 <- lmer(ihs(bio)~Microsite * gams + (1|Year), data=mean.year)
shapiro.test(resid(m1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m1)
## W = 0.9789, p-value = 0.8234
anova(m1)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                 Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Microsite       1.6259  1.6259     1 23.000  11.782  0.002272 ** 
## gams           18.1519 18.1519     1 23.367 131.533 4.448e-11 ***
## Microsite:gams  1.0755  1.0755     1 23.000   7.794  0.010365 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m1) ## R squared values 
##      R2m      R2c 
## 0.549095 0.945669
m2 <- lmer(ihs(abd)~Microsite * gams + (1|Year), data=mean.year)
shapiro.test(resid(m2))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2)
## W = 0.97267, p-value = 0.6537
car::Anova(m2, type=2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ihs(abd)
##                  Chisq Df Pr(>Chisq)    
## Microsite       0.1875  1     0.6650    
## gams           63.2057  1  1.862e-15 ***
## Microsite:gams  0.1083  1     0.7421    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values 
##       R2m       R2c 
## 0.6343625 0.8521801
m2.nat <- lmer(ihs(native)~Microsite * poly(gams,2) + (1|Year), data=mean.status2)
shapiro.test(resid(m2.nat))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2.nat)
## W = 0.98084, p-value = 0.8701
anova(m2.nat)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                         Sum Sq Mean Sq NumDF DenDF F.value  Pr(>F)  
## Microsite               1.2604  1.2604     1    22  1.1237 0.30063  
## poly(gams, 2)           8.0514  4.0257     2    22  3.5890 0.04477 *
## Microsite:poly(gams, 2) 2.5294  1.2647     2    22  1.1275 0.34184  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values 
##       R2m       R2c 
## 0.6343625 0.8521801
m2.exo <- lmer(ihs(non.native)~Microsite * gams+ (1|Year), data=mean.status2)
shapiro.test(resid(m2.exo))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2.exo)
## W = 0.90696, p-value = 0.01672
anova(m2.exo)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Microsite       0.247   0.247     1 23.000   0.183    0.6727    
## gams           74.389  74.389     1 23.965  55.056 1.175e-07 ***
## Microsite:gams  0.175   0.175     1 23.000   0.129    0.7223    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values 
##       R2m       R2c 
## 0.6343625 0.8521801
m2.rat <- lmer(ihs(native.ratio)~Microsite *gams+ (1|Year), data=mean.status2)
shapiro.test(resid(m2.rat))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2.rat)
## W = 0.89603, p-value = 0.009249
anova(m2.rat)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                 Sum Sq Mean Sq NumDF   DenDF F.value  Pr(>F)  
## Microsite      0.03367 0.03367     1 23.0000  0.3040 0.58669  
## gams           0.47126 0.47126     1  9.7629  4.2555 0.06674 .
## Microsite:gams 0.02598 0.02598     1 23.0000  0.2346 0.63271  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values
##       R2m       R2c 
## 0.6343625 0.8521801
m3 <- lmer(rich~Microsite * poly(gams,2) + (1|Year), data=mean.year)
shapiro.test(resid(m3))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m3)
## W = 0.93835, p-value = 0.1004
car::Anova(m3, type=2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: rich
##                           Chisq Df Pr(>Chisq)   
## Microsite                0.6282  1   0.428012   
## poly(gams, 2)           10.9223  2   0.004249 **
## Microsite:poly(gams, 2)  2.5197  2   0.283699   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m3) ## R squared values 
##       R2m       R2c 
## 0.3445392 0.3495475
m4 <- lmer(mpd~Microsite * gams + (1|Year), data=mpd.gams)
shapiro.test(resid(m4))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m4)
## W = 0.98753, p-value = 0.9827
anova(m4)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                 Sum Sq Mean Sq NumDF DenDF F.value   Pr(>F)   
## Microsite      0.68094 0.68094     1    22  9.2604 0.005967 **
## gams           0.27407 0.27407     1    22  3.7272 0.066522 . 
## Microsite:gams 0.55289 0.55289     1    22  7.5189 0.011898 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m4) ## R squared values 
##       R2m       R2c 
## 0.3852134 0.3852134
### Biomass
p1 <- ggplot(mean.year) + geom_jitter(aes(x=gams, y=ihs(bio)), color="#56B4E9", size=3, data=subset(mean.year, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=ihs(bio)), color="#E69F00", size=3, data=subset(mean.year, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,3),aes(x=gams, y=ihs(bio)), color="#E69F00", fill="#E69F0080", data=subset(mean.year, Microsite=="open"), lwd=1.4) +  stat_smooth(method="lm", formula= y~poly(x,3),aes(x=gams, y=ihs(bio)), color="#56B4E9", fill="#56B4E970", data=subset(mean.year, Microsite=="shrub"), lwd=1.4)  + ylab("Annnual plant biomass")+ xlab("")  #annotate("text", x=-4.5,y=4, label="a", size=8)+coord_cartesian(ylim=c(0, 4))

### abundance
p2 <- ggplot(mean.year) + geom_jitter(aes(x=gams, y=ihs(abd)), color="#56B4E9", size=3, data=subset(mean.year, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=ihs(abd)), color="#E69F00", size=3, data=subset(mean.year, Microsite=="open"))  + theme_Publication()+  stat_smooth(method="lm", formula= ihs(y)~x,aes(x=gams, y=abd), color="black", fill="Grey50", data=mean.year, lwd=1.4)+ ylab("Annual plant abundance")+ xlab("")+ annotate("text", x=-4.5,y=6, label="b", size=8)

### richness
p3 <- ggplot(mean.year) + geom_jitter(aes(x=gams, y=rich), color="#56B4E9", size=3, data=subset(mean.year, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=rich), color="#E69F00", size=3, data=subset(mean.year, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=rich), color="#E69F00", fill="#E69F0080", data=subset(mean.year, Microsite=="open"), lwd=1.4) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=rich), color="#56B4E9", fill="#56B4E970", data=subset(mean.year, Microsite=="shrub"), lwd=1.4)  + ylab("Species richness") + xlab("Aridity Gradient")+ annotate("text", x=-4.5,y=4, label="c", size=8)

### phylogenetic diveristy
p4 <- ggplot(mpd.gams)+ geom_jitter(aes(x=gams, y=mpd), color="#56B4E9", size=3, data=subset(mpd.gams, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=mpd), color="#E69F00", size=3, data=subset(mpd.gams, Microsite=="open"))  +  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="shrub"), color="#56B4E9",fill="#56B4E980", lwd=1.4)+  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="open"), color="#E69F00",fill="#E69F0080", lwd=1.4)+ ylab("Phylogenetic Community Dissimilarity") + xlab("Aridity Gradient")+ theme_Publication()+ annotate("text", x=-4.5,y=1.8, label="d", size=8)

grid.arrange(p1,p2,p3,p4)

### Nutrients 
nutrient.mean <- nutrients.climate %>% group_by(Site, arid.gradient, microsite) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean,gams=rep(gams.data[gams.data$Year==2016,"gams"],each=2))

## potassium
m5 <- lm(potassium ~ poly(gams,2)* microsite, data=nutrient.mean)
summary(m5)
## 
## Call:
## lm(formula = potassium ~ poly(gams, 2) * microsite, data = nutrient.mean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.224 -15.052  -5.456   7.166  51.901 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     141.46      12.13  11.662 2.67e-06 ***
## poly(gams, 2)1                 -217.98      45.39  -4.803 0.001351 ** 
## poly(gams, 2)2                  150.46      45.39   3.315 0.010617 *  
## micrositeshrub                   90.20      17.15   5.258 0.000766 ***
## poly(gams, 2)1:micrositeshrub  -256.71      64.19  -4.000 0.003952 ** 
## poly(gams, 2)2:micrositeshrub   149.36      64.19   2.327 0.048389 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.09 on 8 degrees of freedom
## Multiple R-squared:  0.9641, Adjusted R-squared:  0.9416 
## F-statistic: 42.95 on 5 and 8 DF,  p-value: 1.438e-05
anova(m5)
## Analysis of Variance Table
## 
## Response: potassium
##                         Df Sum Sq Mean Sq F value    Pr(>F)    
## poly(gams, 2)            2 170636   85318  82.835 4.502e-06 ***
## microsite                1  28476   28476  27.648 0.0007663 ***
## poly(gams, 2):microsite  2  22053   11026  10.706 0.0054741 ** 
## Residuals                8   8240    1030                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m5) ## R squared values 
##       R2m       R2c 
## 0.9429147 0.9429147
p5 <- ggplot(nutrient.mean) + geom_jitter(aes(x=gams, y=potassium), color="#56B4E9", size=3, data=subset(nutrient.mean, microsite=="shrub"))+ geom_jitter(aes(x=gams, y=potassium), color="#E69F00", size=3, data=subset(nutrient.mean, microsite=="open"))+  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=potassium), color="#56B4E9", fill="#56B4E980", data=subset(nutrient.mean, microsite=="shrub"), lwd=1.4) +  
  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=potassium), color="#E69F00", fill="#E69F0080",data=subset(nutrient.mean, microsite=="open"), lwd=1.4)+ ylab("Nitrogen (ppm)") + xlab("")+ theme_Publication()+ annotate("text", x=-4.5,y=32, label="a", size=8)+coord_cartesian(ylim=c(0, 32))



se <- function(x) sd(na.omit(x))/sqrt(length(na.omit(x)))

smc <- census %>% group_by(Year, Site, Microsite, Census) %>% summarize(swc.avg=mean(swc), swc.se=se(swc))
smc.arid <- merge(smc, gams.data, by=c("Site","Year"))


## early swc
smc.early <- subset(smc.arid, Census=="emergence")


m6 <- lmer(swc.avg~Microsite * poly(gams,2) + (1|Year), data=smc.early)
shapiro.test(resid(m6))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m6)
## W = 0.91432, p-value = 0.02516
anova(m6)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                          Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Microsite                  1.53    1.53     1 21.000   0.084    0.7752    
## poly(gams, 2)           1449.15  724.58     2 21.343  39.560 6.617e-08 ***
## Microsite:poly(gams, 2)    3.69    1.85     2 21.000   0.101    0.9046    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m6) ## R squared values 
##       R2m       R2c 
## 0.5464095 0.8953768
### soil moisture at emergence
p6 <- ggplot(smc.early) + geom_jitter(aes(x=gams, y=swc.avg), color="#56B4E9", size=3, data=subset(smc.early, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=swc.avg), color="#E69F00", size=3, data=subset(smc.early, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=swc.avg), data=smc.early, lwd=1.4, color="black") + ylab("Soil moisture at emergence (%)") + xlab("Aridity Gradient")+ annotate("text", x=-4.5,y=30, label="c", size=8)


hobo.means <- HOBO.data %>% group_by(season, Microsite, Site) %>% summarize(temp.avg=mean(Temp),temp.var=var(Temp),temp.sd=sd(Temp),temp.se=se(Temp),rh.avg=mean(RH, na.rm=T),rh.var=var(RH, na.rm=T),rh.se=se(RH),frost=(sum(na.omit(Temp)<0)/length(na.omit(Temp))*100),heat=(sum(na.omit(Temp)>30)/length(na.omit(Temp))*100))
hobo.means[,"Year"] <- c(rep(2016,14),rep(2017,14))

hobo.arid <- merge(hobo.means, gams.data, by=c("Site","Year"))
# hobo.arid[,"CV"] <-hobo.arid[,"temp.sd"] /hobo.arid[,"CV"] 

## temperature variability
m7 <- lmer(temp.var~Microsite * poly(gams,2) + (1|Year), data=hobo.arid)
shapiro.test(resid(m7))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m7)
## W = 0.96514, p-value = 0.4578
anova(m7)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                         Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Microsite               4406.2  4406.2     1 21.000 21.4200 0.0001448 ***
## poly(gams, 2)           4126.7  2063.3     2 21.483 10.0306 0.0008384 ***
## Microsite:poly(gams, 2)   56.2    28.1     2 21.000  0.1366 0.8730834    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m7) ## R squared values 
##       R2m       R2c 
## 0.4521376 0.7930449
p7 <- ggplot(hobo.arid) + geom_jitter(aes(x=gams, y=temp.var), color="#56B4E9", size=3, data=subset(hobo.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=temp.var), color="#E69F00", size=3, data=subset(hobo.arid, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#E69F00", fill="#E69F0080", data=subset(hobo.arid, Microsite=="open"), lwd=1.4) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#56B4E9", fill="#56B4E970", data=subset(hobo.arid, Microsite=="shrub"), lwd=1.4)  + ylab("Temperature variation") + xlab("Aridity Gradient")+ annotate("text", x=-4.5,y=125, label="d", size=8) 

shrubz <- read.csv("Data/ERG.shrub.csv")
shrub.mean <- shrubz %>% group_by(Site, Microsite) %>% summarize(compact=mean(Compaction))

shrub.clim <- data.frame(shrub.mean, gams=rep(gams.data[gams.data$Year==2016,"gams"],each=2))
#shrub.clim <- subset(shrub.clim, Year==2016)

m8 <- lm(log(compact)  ~  Microsite * gams, data=shrub.clim)
summary(m8)
## 
## Call:
## lm(formula = log(compact) ~ Microsite * gams, data = shrub.clim)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37348 -0.12596  0.04006  0.12843  0.33669 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         -0.138777   0.804662  -0.172   0.8665  
## Micrositeshrub      -3.101869   1.137963  -2.726   0.0213 *
## gams                 0.007791   0.010038   0.776   0.4556  
## Micrositeshrub:gams  0.036028   0.014196   2.538   0.0295 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2383 on 10 degrees of freedom
## Multiple R-squared:  0.6967, Adjusted R-squared:  0.6058 
## F-statistic: 7.658 on 3 and 10 DF,  p-value: 0.006
anova(m8)
## Analysis of Variance Table
## 
## Response: log(compact)
##                Df  Sum Sq Mean Sq F value   Pr(>F)   
## Microsite       1 0.18836 0.18836  3.3176 0.098548 . 
## gams            1 0.75039 0.75039 13.2168 0.004571 **
## Microsite:gams  1 0.36569 0.36569  6.4409 0.029469 * 
## Residuals      10 0.56776 0.05678                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m8) ## R squared values 
##       R2m       R2c 
## 0.6386405 0.6386405
p8 <- ggplot(shrub.clim) + geom_jitter(aes(x=gams, y=compact), color="#56B4E9", size=3, data=subset(shrub.clim, Microsite=="shrub"))+ geom_jitter(aes(x=gams, y=compact), color="#E69F00", size=3, data=subset(shrub.clim, Microsite=="open"))+  stat_smooth(method="lm", formula= y~exp(x),aes(x=gams, y=compact), color="#56B4E9", fill="#56B4E980", data=subset(shrub.clim, Microsite=="shrub"), lwd=1.4) +  
  stat_smooth(method="lm", formula= y~exp(x),aes(x=gams, y=compact), color="#E69F00", fill="#E69F0080",data=subset(shrub.clim, Microsite=="open"), lwd=1.4)+ ylab("Soil compaction") + xlab("")+ theme_Publication() + annotate("text", x=-4.4,y=2.5, label="b", size=8) +coord_cartesian(ylim=c(0.5, 2.5))

grid.arrange(p5,p8,p6,p7)

Microsite comparisons

nut <- read.csv("Data//ERG.soilnutrients.csv")

fit1 <- aovp(N ~ microsite, data=nut)
## [1] "Settings:  unique SS "
summary(fit1)
## Component 1 :
##             Df R Sum Sq R Mean Sq Iter Pr(Prob)   
## microsite    1   1031.1    1031.1 5000   0.0014 **
## Residuals   68   6582.7      96.8                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t.test(N ~ microsite, data=nut)
## 
##  Welch Two Sample t-test
## 
## data:  N by microsite
## t = -3.2637, df = 36.403, p-value = 0.002396
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -12.444151  -2.907849
## sample estimates:
##  mean in group open mean in group shrub 
##            1.594571            9.270571
fit2 <- aovp(P ~ microsite, data=nut)
## [1] "Settings:  unique SS "
t.test(P ~ microsite, data=nut)
## 
##  Welch Two Sample t-test
## 
## data:  P by microsite
## t = -2.2393, df = 53.804, p-value = 0.02929
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.094126 -0.391588
## sample estimates:
##  mean in group open mean in group shrub 
##            6.365714           10.108571
summary(fit2)
## Component 1 :
##             Df R Sum Sq R Mean Sq Iter Pr(Prob)  
## microsite    1    245.2   245.157 5000   0.0106 *
## Residuals   68   3324.4    48.889                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit3 <- aovp(K ~ microsite, data=nut)
## [1] "Settings:  unique SS "
t.test(K ~ microsite, data=nut)
## 
##  Welch Two Sample t-test
## 
## data:  K by microsite
## t = -2.8179, df = 55.964, p-value = 0.006669
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -154.32488  -26.07512
## sample estimates:
##  mean in group open mean in group shrub 
##            141.4571            231.6571
summary(fit3)
## Component 1 :
##             Df R Sum Sq R Mean Sq Iter Pr(Prob)   
## microsite    1   142381    142381 5000    0.003 **
## Residuals   68  1219331     17931                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
compact <- read.csv("Data//ERG.shrub.csv")

fit4 <- aovp(Compaction ~ Microsite, data=compact)
## [1] "Settings:  unique SS "
t.test(Compaction ~ Microsite, data=compact)
## 
##  Welch Two Sample t-test
## 
## data:  Compaction by Microsite
## t = 2.6675, df = 411.15, p-value = 0.007943
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.06357828 0.41975506
## sample estimates:
##  mean in group open mean in group shrub 
##            1.653571            1.411905
summary(fit4)
## Component 1 :
##              Df R Sum Sq R Mean Sq Iter Pr(Prob)   
## Microsite     1     6.13    6.1323 5000   0.0044 **
## Residuals   418   360.23    0.8618                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
soilmoist <- read.csv("Data//ERG.phytometer.census.csv")

fit5 <- aovp(swc ~ Microsite, data=subset(soilmoist, Census=="emergence"))
## [1] "Settings:  unique SS "
summary(fit5)
## Component 1 :
##              Df R Sum Sq R Mean Sq Iter Pr(Prob)
## Microsite     1       46    45.967  179   0.3631
## Residuals   838    62917    75.080
t.test(swc ~ Microsite, data=subset(soilmoist, Census=="emergence"))
## 
##  Welch Two Sample t-test
## 
## data:  swc by Microsite
## t = 0.78246, df = 827.4, p-value = 0.4342
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.7057865  1.6415008
## sample estimates:
##  mean in group open mean in group shrub 
##            11.16500            10.69714
micro <- read.csv("Data//ERG.logger.data.csv")

micro.avg <- micro %>% group_by(Year, Site, Microsite, Rep) %>% summarize(temp=mean(Temp), rh=mean(RH))

fit6 <- aovp(temp ~ Microsite, data=micro.avg)
## [1] "Settings:  unique SS "
summary(fit6)
## Component 1 :
##              Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## Microsite     1    77.61    77.608 5000 < 2.2e-16 ***
## Residuals   111   903.55     8.140                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t.test(temp ~ Microsite, data=micro.avg)
## 
##  Welch Two Sample t-test
## 
## data:  temp by Microsite
## t = 3.0867, df = 110.64, p-value = 0.002558
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.5933932 2.7216712
## sample estimates:
##  mean in group open mean in group shrub 
##           11.101469            9.443937
fit7 <- aovp(rh ~ Microsite, data=micro.avg)
## [1] "Settings:  unique SS "
summary(fit7)
## Component 1 :
##              Df R Sum Sq R Mean Sq Iter Pr(Prob)
## Microsite     1      375    375.35  114   0.4737
## Residuals   108    34999    324.07
t.test(rh ~ Microsite, data=micro.avg)
## 
##  Welch Two Sample t-test
## 
## data:  rh by Microsite
## t = -1.0736, df = 105.97, p-value = 0.2854
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -10.523642   3.129837
## sample estimates:
##  mean in group open mean in group shrub 
##            55.04511            58.74201

Plot level GLMM

m1 <- glmer.nb(Abundance ~ Microsite * gams + (1|Year), data=community.arid)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
car::Anova(m1, type=2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Abundance
##                    Chisq Df Pr(>Chisq)    
## Microsite        16.3613  1  5.234e-05 ***
## gams           1136.3911  1  < 2.2e-16 ***
## Microsite:gams    1.0278  1     0.3107    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef(m1)
## $Year
##      (Intercept) Micrositeshrub       gams Micrositeshrub:gams
## 2016    13.13187      0.6544052 -0.1298036        -0.005782823
## 2017    11.53434      0.6544052 -0.1298036        -0.005782823
## 
## attr(,"class")
## [1] "coef.mer"
ggplot(community.arid) + geom_jitter(aes(x=gams, y=Abundance), width = 0.5, color="#56B4E9", size=3, data=subset(community.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=Abundance),width = 0.5, color="#E69F00", size=3, data=subset(community.arid, Microsite=="open"))  + theme_Publication() + stat_function(fun = function(x) exp(12.332+ x*-0.1298), lwd=2, lty=1, col="Grey40")+  ylab("annual plant abundance") + xlab("Gams index of continentality")
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

source("rsquared.r")

m2 <- glmer.nb(Richness ~ Microsite * poly(gams,2) + (1|Year), data=community.arid)
## Warning in theta.ml(Y, mu, weights = object@resp$weights, limit = limit, :
## iteration limit reached
car::Anova(m2, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Richness
##                            Chisq Df Pr(>Chisq)    
## (Intercept)             140.6109  1  < 2.2e-16 ***
## Microsite                 8.7643  1   0.003072 ** 
## poly(gams, 2)           969.9903  2  < 2.2e-16 ***
## Microsite:poly(gams, 2) 374.2555  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# rsquared.glmm(m2)

newdf <- expand.grid(gams=seq(45,90, by=0.05), Microsite=c("shrub","open"))
newdf$fit = predict(m2, type = "response", newdata = newdf)
## Warning in pred + base::drop(as(newRE$b %*% newRE$Zt, "matrix")): longer
## object length is not a multiple of shorter object length
## richness
p1 <- ggplot(community.arid) + geom_jitter(aes(x=gams, y=Richness), color="#181818", size=3, alpha=0.4, width = 0.75,shape = 1, data=subset(community.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=Richness), color="#181818", size=2,alpha=0.4, data=subset(community.arid, Microsite=="open"),shape = 2, width = 0.75)  + theme_Publication() +  ylab("species richness") + xlab("Gams index of continentality")+ 
  geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~poly(x,2),aes(x=gams, y=Richness), color="#181818", fill="#80808080" , data=subset(community.arid, Microsite=="shrub"), lwd=2) + 
  geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~poly(x,2),aes(x=gams, y=Richness), color="#181818", fill="#80808080", data=subset(community.arid, Microsite=="open"), lwd=2, lty=2)+  ylab("Species richness") + xlab("") + annotate("text", x=50,y=6, label="a", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
#+ geom_smooth(data=subset(newdf,Microsite=="shrub"), aes(x=gams, y=fit), color="#181818",lwd=2)+ geom_smooth(data=subset(newdf,Microsite=="open"), aes(x=gams, y=fit),color="#505050",lwd=2)  

m3<- lmer(ihs(Biomass) ~ Microsite * poly(gams,2) + (1|Year), data=community.arid)
anova(m3)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                         Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Microsite                64.53  64.527     1 833.00  210.49 < 2.2e-16 ***
## poly(gams, 2)           511.61 255.805     2 833.44  834.46 < 2.2e-16 ***
## Microsite:poly(gams, 2)  32.08  16.041     2 833.00   52.33 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m3) ## R squared values
##       R2m       R2c 
## 0.5031296 0.8793048
## biomass
p2 <- ggplot(community.arid) + geom_jitter(aes(x=gams, y=ihs(Biomass)), color="#181818", size=3, width = 0.75,shape = 1, alpha=0.5, data=subset(community.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=ihs(Biomass)), color="#181818", size=2,shape = 2,alpha=0.5, width = 0.75, data=subset(community.arid, Microsite=="open"))  + theme_Publication() + 
  geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=ihs(Biomass)), color="#181818", fill="#80808080", data=subset(community.arid, Microsite=="shrub"), lwd=2) + 
  geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=ihs(Biomass)), color="#181818", fill="#80808080", data=subset(community.arid, Microsite=="open"), lwd=2, lty=2)+  ylab("Annual biomass") + xlab("")+ annotate("text", x=50,y=5, label="b", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
m4<- lmer(mpd~ Microsite * gams + (1|Year), data=mpd.gams)
anova(m4)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                 Sum Sq Mean Sq NumDF DenDF F.value   Pr(>F)   
## Microsite      0.68094 0.68094     1    22  9.2604 0.005967 **
## gams           0.27407 0.27407     1    22  3.7272 0.066522 . 
## Microsite:gams 0.55289 0.55289     1    22  7.5189 0.011898 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m4) ## R squared values
##       R2m       R2c 
## 0.3852134 0.3852134
newdf <- expand.grid(gams=seq(75,90, by=0.05), Microsite=c("shrub","open"))
newdf$fit = predict(m4, type = "response", newdata = newdf)

## phylogenetic diversity
p3 <- ggplot(mpd.gams)+ geom_jitter(aes(x=gams, y=mpd), color="#181818", size=3,shape=1, alpha=0.4, data=subset(mpd.gams, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=mpd), color="#181818",alpha=0.4, size=2,shape=2, data=subset(mpd.gams, Microsite=="open"))  +  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="shrub"), color="#181818",fill="#80808080", lwd=2)+  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="open"), color="#181818",fill="#80808080", lwd=2, lty=2)+ ylab("Phylogenetic community dissimilarity") + xlab("Gams index of continentality")+ theme_Publication()+ annotate("text", x=50,y=1.4, label="c", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## native vs non-native
status.gams <- merge(status.data, gams.data, by=c("Year","Site"))

m5<- glmer.nb(native~ Microsite * scale(gams) + (1|Year), data=status.gams) ## need to scale data
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00214781 (tol =
## 0.001, component 1)
car::Anova(m5, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: native
##                         Chisq Df Pr(>Chisq)    
## (Intercept)           565.724  1  < 2.2e-16 ***
## Microsite              23.122  1  1.521e-06 ***
## scale(gams)            67.393  1  2.224e-16 ***
## Microsite:scale(gams)  35.708  1  2.292e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# rsquared.glmm(m5)


m6<- glmer.nb(non.native~ Microsite * scale(gams) + (1|Year), data=status.gams)   ## need to scale data
car::Anova(m6, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: non.native
##                          Chisq Df Pr(>Chisq)    
## (Intercept)            12.5323  1   0.000400 ***
## Microsite              13.5445  1   0.000233 ***
## scale(gams)           353.8669  1  < 2.2e-16 ***
## Microsite:scale(gams)   1.4368  1   0.230664    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# rsquared.glmm(m5)


status.plot <- gather(status.gams, origin, abundance, 5:6)

p4 <- ggplot(subset(status.plot, abundance>0)) + geom_jitter(aes(x=gams, y=abundance), color="#181818", size=3, data=subset(status.plot, Microsite=="shrub"), width = 1, shape=1, alpha=0.3) +geom_jitter(aes(x=gams, y=abundance), width = 1, color="#181818", size=2,alpha=0.3,shape=2, data=subset(status.plot, Microsite=="open"))  + theme_Publication() + 
  geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance), color="#505050", se=F, data=subset(status.plot, Microsite=="open" & origin=="native"), lwd=2, lty=2)+
  geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance),se=F, color="#181818",  data=subset(status.plot, Microsite=="open" & origin=="non.native"), lwd=2, lty=2)+
 geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance),se=F, color="#505050", data=subset(status.plot, Microsite=="shrub" & origin=="native"), lwd=2)+
  geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance),se=F, color="#181818",  data=subset(status.plot, Microsite=="shrub" & origin=="non.native"), lwd=2) +  ylab("Annual plant abundance") + xlab("Gams index of continentality")+ annotate("text", x=50,y=310, label="d", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
grid.arrange(p1,p2,p3,p4)
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

### Nutrients 
nutrient.mean <- nutrients.climate %>% group_by(Site, arid.gradient, microsite) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean,gams=rep(gams.data[gams.data$Year==2016,"gams"],each=2))

## nitrogen
m5 <- lm(nitrogen ~ poly(gams,2)* microsite, data=nutrient.mean)
summary(m5)
## 
## Call:
## lm(formula = nitrogen ~ poly(gams, 2) * microsite, data = nutrient.mean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2918 -0.9251 -0.2773  0.8012  5.6374 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1.595      1.038   1.536 0.163115    
## poly(gams, 2)1                   3.927      3.885   1.011 0.341676    
## poly(gams, 2)2                   3.847      3.885   0.990 0.350969    
## micrositeshrub                   7.676      1.468   5.228 0.000795 ***
## poly(gams, 2)1:micrositeshrub   22.258      5.494   4.052 0.003676 ** 
## poly(gams, 2)2:micrositeshrub   17.833      5.494   3.246 0.011771 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.747 on 8 degrees of freedom
## Multiple R-squared:  0.9298, Adjusted R-squared:  0.8859 
## F-statistic: 21.18 on 5 and 8 DF,  p-value: 0.0002012
anova(m5)
## Analysis of Variance Table
## 
## Response: nitrogen
##                         Df Sum Sq Mean Sq F value    Pr(>F)    
## poly(gams, 2)            2 389.59 194.794  25.817 0.0003239 ***
## microsite                1 206.22 206.223  27.332 0.0007948 ***
## poly(gams, 2):microsite  2 203.35 101.677  13.476 0.0027446 ** 
## Residuals                8  60.36   7.545                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m5) ## R squared values 
##       R2m       R2c 
## 0.8906817 0.8906817
p5 <- ggplot(nutrient.mean) + geom_jitter(aes(x=gams, y=nitrogen), color="#181818", size=3, shape=16,alpha=0.5, data=subset(nutrient.mean, microsite=="shrub"))+ geom_jitter(aes(x=gams, y=nitrogen), color="#181818", size=3,shape=17, alpha=0.4, data=subset(nutrient.mean, microsite=="open"))+  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=nitrogen), color="#181818", fill="#80808080", data=subset(nutrient.mean, microsite=="shrub"), lwd=2) +  
  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=nitrogen), color="#181818", fill="#80808080",data=subset(nutrient.mean, microsite=="open"), lwd=2, lty=2)+ ylab("Nitrogen (ppm)") + xlab("")+ theme_Publication()+ coord_cartesian(ylim=c(-2, 26)) + xlab("")+ annotate("text", x=63,y=25, label="a", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
smc.early <- subset(smc, Census=="emergence")
smc.early <- merge(smc.early, gams.data, by=c("Year","Site"))

## nitrogen
m6 <- lmer(log(swc.avg) ~ poly(gams,2)* Microsite + (1|Year), data=smc.early)
summary(m6)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: log(swc.avg) ~ poly(gams, 2) * Microsite + (1 | Year)
##    Data: smc.early
## 
## REML criterion at convergence: 16.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.75790 -0.67057  0.05989  0.64573  1.52776 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Year     (Intercept) 0.32730  0.5721  
##  Residual             0.09642  0.3105  
## Number of obs: 28, groups:  Year, 2
## 
## Fixed effects:
##                                Estimate Std. Error        df t value
## (Intercept)                    2.183930   0.412962  1.000000   5.288
## poly(gams, 2)1                -4.095881   0.529751 21.483000  -7.732
## poly(gams, 2)2                 0.009355   0.444206 21.040000   0.021
## Micrositeshrub                -0.021937   0.117363 21.000000  -0.187
## poly(gams, 2)1:Micrositeshrub  0.108296   0.621024 21.000000   0.174
## poly(gams, 2)2:Micrositeshrub -0.128193   0.621024 21.000000  -0.206
##                               Pr(>|t|)    
## (Intercept)                      0.119    
## poly(gams, 2)1                1.22e-07 ***
## poly(gams, 2)2                   0.983    
## Micrositeshrub                   0.854    
## poly(gams, 2)1:Micrositeshrub    0.863    
## poly(gams, 2)2:Micrositeshrub    0.838    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(,2)1 pl(,2)2 Mcrsts p(,2)1:
## ply(gms,2)1  0.000                               
## ply(gms,2)2  0.000  0.084                        
## Microstshrb -0.142  0.000   0.000                
## ply(g,2)1:M  0.000 -0.586   0.000   0.000        
## ply(g,2)2:M  0.000  0.000  -0.699   0.000  0.000
anova(m6)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                         Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## poly(gams, 2)           8.6913  4.3456     2 21.339  45.071 2.183e-08 ***
## Microsite               0.0034  0.0034     1 21.000   0.035    0.8535    
## poly(gams, 2):Microsite 0.0070  0.0035     2 21.000   0.037    0.9642    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m6) ## R squared values 
##       R2m       R2c 
## 0.5883174 0.9063213
p6 <- ggplot(smc.early) + geom_jitter(aes(x=gams, y=swc.avg), color="#181818", size=3, shape=16,alpha=0.5, data=subset(smc.early, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=swc.avg), color="#181818", size=3,shape=17,alpha=0.5, data=subset(smc.early, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=swc.avg), data=smc.early, lwd=2, color="black") + ylab("Soil moisture at emergence (%)")  + xlab("")+ annotate("text", x=50,y=35, label="b", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## temperature variability
m7 <- lmer(temp.var~Microsite * poly(gams,2) + (1|Year), data=hobo.arid)
shapiro.test(resid(m7))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m7)
## W = 0.96514, p-value = 0.4578
anova(m7)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                         Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
## Microsite               4406.2  4406.2     1 21.000 21.4200 0.0001448 ***
## poly(gams, 2)           4126.7  2063.3     2 21.483 10.0306 0.0008384 ***
## Microsite:poly(gams, 2)   56.2    28.1     2 21.000  0.1366 0.8730834    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m7) ## R squared values 
##       R2m       R2c 
## 0.4521376 0.7930449
p7 <- ggplot(hobo.arid) + geom_jitter(aes(x=gams, y=temp.var), color="#181818", size=3, shape=16, alpha=0.5, data=subset(hobo.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=temp.var), color="#181818", size=3, shape=17, alpha=0.5, data=subset(hobo.arid, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#181818", fill="#80808080", data=subset(hobo.arid, Microsite=="open"), lwd=2) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#181818", fill="#80808080", data=subset(hobo.arid, Microsite=="shrub"), lwd=2, lty=2)  + ylab("Temperature variation") + xlab("Gams index of continentality")+ annotate("text", x=50,y=125, label="c", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## soil compaction
m8 <- lm(log(compact)  ~  Microsite * gams, data=shrub.clim)
summary(m8)
## 
## Call:
## lm(formula = log(compact) ~ Microsite * gams, data = shrub.clim)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37348 -0.12596  0.04006  0.12843  0.33669 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         -0.138777   0.804662  -0.172   0.8665  
## Micrositeshrub      -3.101869   1.137963  -2.726   0.0213 *
## gams                 0.007791   0.010038   0.776   0.4556  
## Micrositeshrub:gams  0.036028   0.014196   2.538   0.0295 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2383 on 10 degrees of freedom
## Multiple R-squared:  0.6967, Adjusted R-squared:  0.6058 
## F-statistic: 7.658 on 3 and 10 DF,  p-value: 0.006
anova(m8)
## Analysis of Variance Table
## 
## Response: log(compact)
##                Df  Sum Sq Mean Sq F value   Pr(>F)   
## Microsite       1 0.18836 0.18836  3.3176 0.098548 . 
## gams            1 0.75039 0.75039 13.2168 0.004571 **
## Microsite:gams  1 0.36569 0.36569  6.4409 0.029469 * 
## Residuals      10 0.56776 0.05678                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m8) ## R squared values 
##       R2m       R2c 
## 0.6386405 0.6386405
p8 <- ggplot(shrub.clim) + geom_jitter(aes(x=gams, y=log(compact)), color="#181818", size=3, shape=16, alpha=0.5, data=subset(shrub.clim, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=log(compact)), color="#181818", size=3, shape=17, alpha=0.5, data=subset(shrub.clim, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~x, aes(x=gams, y=log(compact)), color="#181818", fill="#80808080", data=subset(shrub.clim, Microsite=="open"), lwd=2) +  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=log(compact)), color="#181818", fill="#80808080", data=subset(shrub.clim, Microsite=="shrub"), lwd=2, lty=2)  + ylab("Soil compaction") + xlab("Gams index of continentality")+ annotate("text", x=63,y=1.1, label="d", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
grid.arrange(p5,p6,p7,p8)
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

Phytometer with GAMS

## treatment plots
census.long <- gather(census.arid, species, biomass, 12:14)
census.long[,"species"] <- ifelse(census.long[,"species"] =="Phacelia.biomass", "P. tanacetifolia",
                                  ifelse(census.long[,"species"]=="Plantago.biomass","P. insularis","S. columbariae"))
census.long[,"species"] <- as.factor(census.long$species)

## microsite average
census.micro <- census.long %>% filter(biomass>0) %>%  group_by(Microsite, species) %>%  summarize(Biomass=mean(biomass), error=se(biomass))

p9 <- ggplot(census.micro, aes(x=species, y=Biomass, fill=Microsite)) +   geom_bar(position=position_dodge(), stat="identity", colour="black") + scale_fill_brewer(palette="Greys")+  geom_errorbar(aes(ymin=Biomass-error, ymax=Biomass+error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+ theme_Publication() + theme(axis.text.x = element_text(size=16, face="italic"),axis.text.y = element_text(size=16),axis.title = element_text(size=20),legend.text=element_text(size=18))
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## nutrient average
census.nut <- census.long %>% filter(biomass>0) %>%  group_by(Nutrient, species) %>%  summarize(Biomass=mean(biomass), error=se(biomass))

p10 <- ggplot(census.nut, aes(x=species, y=Biomass, fill=Nutrient)) +   geom_bar(position=position_dodge(), stat="identity",colour="black") + scale_fill_brewer(palette="Greys")+  geom_errorbar(aes(ymin=Biomass-error, ymax=Biomass+error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+ theme_Publication()+ theme(axis.text.x = element_text(size=16, face="italic"),axis.text.y = element_text(size=16),axis.title = element_text(size=20),legend.text=element_text(size=18))
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## aridity average

# ggplot(census.long, aes(x=arid.gradient, y=biomass, fill=species, colour=species)) + theme_Publication()+   stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=biomass), color="#009E73", fill="#009E7340", data=subset(census.long, species=="Phacelia.biomass" & biomass>0)) + stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=biomass), color="#505050",fill="#50505040", data=subset(census.long, species=="Plantago.biomass" & biomass>0))+ stat_smooth(method="lm", formula= y~exp(x),aes(x=arid.gradient, y=biomass), color="#181818", fill="#18181840", data=subset(census.long, species=="Salvia.biomass" & biomass>0))

### need to separate zeros using gamma-hurdle model
##http://seananderson.ca/2014/05/18/gamma-hurdle.html


## logistic first
census.long[,"presence"] <- ifelse(census.long$biomass>0,1,0)
census.long <- merge(census.long, gams.data, by=c("Year","Site"))
## Drop panoche 2017 because of cold
census.npan <- census.long %>% filter(Site!="PanocheHills" | Year!="2017")


m1.l <- glmer(presence~Microsite * gams * Nutrient + (1|Year),data=subset(census.npan, species=="P. tanacetifolia"), family=binomial)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
car::Anova(m1.l, type=2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: presence
##                           Chisq Df Pr(>Chisq)    
## Microsite               29.2566  1   6.34e-08 ***
## gams                     8.9859  1   0.002721 ** 
## Nutrient                 0.0060  1   0.938002    
## Microsite:gams           4.4524  1   0.034853 *  
## Microsite:Nutrient       1.9750  1   0.159917    
## gams:Nutrient            0.0019  1   0.965645    
## Microsite:gams:Nutrient  0.0139  1   0.906248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- lmer(log(biomass)~Microsite * poly(gams,2) * Nutrient + (1|Year),data=subset(census.long, species=="P. tanacetifolia" & presence==1))
shapiro.test(resid(m1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m1)
## W = 0.99225, p-value = 0.2595
anova(m1)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                                   Sum Sq Mean Sq NumDF DenDF F.value
## Microsite                         22.121  22.121     1   221  12.213
## poly(gams, 2)                    157.099  78.549     2   221  43.368
## Nutrient                           5.710   5.710     1   221   3.152
## Microsite:poly(gams, 2)           15.688   7.844     2   221   4.331
## Microsite:Nutrient                 0.078   0.078     1   221   0.043
## poly(gams, 2):Nutrient             2.854   1.427     2   221   0.788
## Microsite:poly(gams, 2):Nutrient   1.406   0.703     2   221   0.388
##                                     Pr(>F)    
## Microsite                        0.0005733 ***
## poly(gams, 2)                    < 2.2e-16 ***
## Nutrient                         0.0771914 .  
## Microsite:poly(gams, 2)          0.0142925 *  
## Microsite:Nutrient               0.8361223    
## poly(gams, 2):Nutrient           0.4561215    
## Microsite:poly(gams, 2):Nutrient 0.6786980    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m1) ## R squared values 
##      R2m      R2c 
## 0.346015 0.346015
m2.l <- glmer(presence~Microsite * poly(gams,2) * Nutrient + (1|Year),data=subset(census.npan, species=="P. insularis"), family=binomial)
car::Anova(m2.l, type=2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: presence
##                                   Chisq Df Pr(>Chisq)   
## Microsite                        1.6804  1   0.194870   
## poly(gams, 2)                    9.6027  2   0.008219 **
## Nutrient                         0.3064  1   0.579874   
## Microsite:poly(gams, 2)          5.3976  2   0.067285 . 
## Microsite:Nutrient               1.9903  1   0.158314   
## poly(gams, 2):Nutrient           3.4679  2   0.176586   
## Microsite:poly(gams, 2):Nutrient 0.6252  2   0.731534   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2 <- lmer(log(biomass)~Microsite * poly(gams,2) * Nutrient + (1|Year), data=subset(census.long, species=="P. insularis" & presence==1))
shapiro.test(resid(m2))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2)
## W = 0.99193, p-value = 0.5364
anova(m2)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                                  Sum Sq Mean Sq NumDF  DenDF F.value
## Microsite                         0.569  0.5693     1 141.23  0.5147
## poly(gams, 2)                    36.153 18.0765     2 137.31 16.3429
## Nutrient                          8.958  8.9576     1 141.27  8.0986
## Microsite:poly(gams, 2)           0.858  0.4289     2 141.51  0.3878
## Microsite:Nutrient                0.110  0.1096     1 141.02  0.0991
## poly(gams, 2):Nutrient            2.373  1.1866     2 141.20  1.0728
## Microsite:poly(gams, 2):Nutrient  0.654  0.3272     2 141.03  0.2958
##                                     Pr(>F)    
## Microsite                          0.47430    
## poly(gams, 2)                    4.297e-07 ***
## Nutrient                           0.00509 ** 
## Microsite:poly(gams, 2)            0.67929    
## Microsite:Nutrient                 0.75338    
## poly(gams, 2):Nutrient             0.34481    
## Microsite:poly(gams, 2):Nutrient   0.74439    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values 
##       R2m       R2c 
## 0.2802814 0.5695942
m3.l <- glmer(presence~Microsite * poly(gams,2) * Nutrient + (1|Year),data=subset(census.npan, species=="S. columbariae"), family=binomial)
car::Anova(m3.l, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: presence
##                                    Chisq Df Pr(>Chisq)    
## (Intercept)                       2.1914  1    0.13878    
## Microsite                         1.6046  1    0.20525    
## poly(gams, 2)                    25.8433  2  2.445e-06 ***
## Nutrient                          0.1181  1    0.73109    
## Microsite:poly(gams, 2)           1.8388  2    0.39875    
## Microsite:Nutrient                3.3815  1    0.06593 .  
## poly(gams, 2):Nutrient            0.9600  2    0.61878    
## Microsite:poly(gams, 2):Nutrient  0.4124  2    0.81369    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 <- lmer(log(biomass)~Microsite * poly(gams,2) * Nutrient + (1|Year), data=subset(census.long, species=="S. columbariae" & biomass>0))
shapiro.test(resid(m3))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m3)
## W = 0.99363, p-value = 0.1532
anova(m3)
## Analysis of Variance Table of type III  with  Satterthwaite 
## approximation for degrees of freedom
##                                  Sum Sq Mean Sq NumDF  DenDF F.value
## Microsite                         0.350  0.3503     1 333.03  0.2070
## poly(gams, 2)                    61.028 30.5139     2 327.98 18.0272
## Nutrient                         18.797 18.7973     1 333.12 11.1052
## Microsite:poly(gams, 2)           1.304  0.6521     2 333.03  0.3853
## Microsite:Nutrient                0.386  0.3859     1 333.05  0.2280
## poly(gams, 2):Nutrient            2.914  1.4569     2 333.05  0.8607
## Microsite:poly(gams, 2):Nutrient  3.356  1.6780     2 333.01  0.9913
##                                     Pr(>F)    
## Microsite                        0.6494522    
## poly(gams, 2)                    3.733e-08 ***
## Nutrient                         0.0009576 ***
## Microsite:poly(gams, 2)          0.6805772    
## Microsite:Nutrient               0.6333374    
## poly(gams, 2):Nutrient           0.4237872    
## Microsite:poly(gams, 2):Nutrient 0.3721737    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m3) ## R squared values 
##       R2m       R2c 
## 0.1506253 0.2337409
## summary presence
mean.occ <- census.npan %>% group_by(gams,species, Microsite) %>% summarize(presence=mean(presence))


p11 <- ggplot(census.npan, aes(x=gams, y=presence, fill=species, color=species))  + theme_Publication()+  
  geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~x,aes(x=gams, y=presence),fill="#0daec2", color="#0daec2",  data=subset(census.npan, species=="P. tanacetifolia" & Microsite=="shrub"), lwd=2)+  
  geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~x,aes(x=gams, y=presence),fill="#0daec2", color="#0daec2",  data=subset(census.npan, species=="P. tanacetifolia"& Microsite=="open"), lwd=2, lty=2) +  geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~poly(x,2),aes(x=gams, y=presence), color="#fb6f70", fill="#fb6f7050",  data=subset(census.npan, species=="P. insularis"), lwd=2)+ geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~poly(x,2),aes(x=gams, y=presence),fill="#b7d73850", color="#b7d738",  data=subset(census.npan, species=="S. columbariae"), lwd=2)+ ylab("Probability of occurrence") + xlab("Gams index of continentality")
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
p12 <- ggplot(census.long, aes(x=gams, y=biomass, fill=species, color=species))  + theme_Publication() +  geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)), color="#fb6f70", fill="#fb6f7050",  data=subset(census.long, species=="P. insularis" & biomass>0), lwd=2, fullrange=TRUE) + geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)),fill="#b7d73850", color="#b7d738",  data=subset(census.long, species=="S. columbariae"& biomass>0), lwd=3)+
  geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)),fill="#0daec2", color="#0daec2",  data=subset(census.long, species=="P. tanacetifolia" & biomass>0 & Microsite=="shrub"), lwd=2, fullrange=TRUE)+geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)),fill="#0daec2", color="#0daec2",  data=subset(census.long, species=="P. tanacetifolia" & biomass>0 & Microsite=="open"), lwd=2,lty=2, fullrange=TRUE)+ ylab("Phytometer biomass") + xlab("Gams index of continentality") 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
grid.arrange(p9,p11,p10, p12)
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## test to see if GAMS related to aridity

arid.gam <- merge(arid.vals, gams.data, by=c("Site"))

cor(arid.gam$Gradient, arid.gam$gams)
## [1] -0.5506378